I am trying to solve equation in the form of y'' + ay' + by + c = 0 (second order differential equation) in python using odeint. As I understood, odeint works only with initial conditions in the form of y(0) = y1, y'(0) = y2. Whereas my conditions are boundary: y'(0) = 0, y'(pi/4) = 0. Is there any way to use odeint with such conditions? Here is the code I have:
from scipy.integrate import odeint
import numpy as np
def lagrange(x, teta):
y = x[0]
dy = x[1]
xdot = [[],[]]
xdot[0] = dy
xdot[1] = (y + dy**2*y*(y**2 + dy**2)**(-1.5) + y*(y**2+dy**2)**(-0.5))/((y**2+dy**2)**(-0.5) - dy**2*(y**2+dy**2)**(-1.5))
return xdot
phi = np.linspace(0,np.pi/4,100)
U = odeint(lagrange, [u1_0, u2_0], phi)
The official documentation states that it solves the initial value problem for stiff or non-stiff systems of first order ode-s. Therefore it cannot solve boundary value problems.
You will most likely need implement the finite difference method, finite element method or shooting method to solve the problem.
I would personally recommend the finite difference method as it essentially consists of adding 3 diagonal matrics (2 of which are tridiagonal), setting up a solution vector and using a linear solver.
If you need some help setting it up, I wrote a diffusion-advection solver not too long ago for matlab which you could port and extend to include the reaction term (should be b*identity matrix).
Depending on your value of a, you may also need to implement artificial diffusion (in case that the peclet number is too large).