I have a system of 3 differential equations (will be obvious from the code I believe) with 3 boundary conditions. I managed to solve it in MATLAB with a loop to change the initial guess bit by bit without terminating the program if it is about to return an error. However, on `scipy`

's `solve_bvp`

, I always get *some* answer, although it is wrong. So I kept changing my guesses (which kept changing the answer) and am giving pretty close numbers to what I have from the actual solution and it's still not working. Is there some other problem with the code perhaps, due to which it's not working? I just edited their documentation's code.

```
import numpy as np
def fun(x, y):
return np.vstack((3.769911184e12*np.exp(-19846/y[1])*(1-y[0]), 0.2056315191*(y[2]-y[1])+6.511664773e14*np.exp(-19846/y[1])*(1-y[0]), 1.696460033*(y[2]-y[1])))
def bc(ya, yb):
return np.array([ya[0], ya[1]-673, yb[2]-200])
x = np.linspace(0, 1, 5)
#y = np.ones((3, x.size))
y = np.array([[1, 1, 1, 1, 1], [670, 670, 670, 670, 670], [670, 670, 670, 670, 670] ])
from scipy.integrate import solve_bvp
sol = solve_bvp(fun, bc, x, y)
```

The actual solution is given below in the figure.

MATLAB Solution to the BVP

## Best Solution

Apparently you need a better initial guess, otherwise the iterative method used by

`solve_bvp`

can create values in`y[1]`

that make the expression`exp(-19846/y[1])`

overflow. When that happens, the algorithm is likely to fail. An overflow in that expression means that some value in`y[1]`

is negative; i.e. the solver is so far out in the weeds that it has little chance of converging to a correct solution. You'll see warnings, and sometimes the function still returns a reasonable solution, but usually it returns garbage when the overflow occurs.You can determine whether or not

`solve_bvp`

has failed to converge by checking`sol.status`

. If it is not 0, something failed.`sol.message`

contains a text message describing the status.I was able to get the Matlab solution by using this to create the initial guess:

Smaller values of

`n`

also work, but when`n`

is too small, an overflow warning can appear.Here's my modified version of your script, followed by the plot that it generates: