# Solve a system of coupled differential equations in Python

I have a system of two coupled differential equations, one is a third-order and the second is second-order. I am looking for a way to solve it in Python.

I would be extremely grateful for any advice on how can I do that or simplify this set of equations that define a boundary value problem :

Pr is just a constant (Prandtl number)

N.B : This question is not related to any previous topic, this set of equation may need a simplification that I don't see. Thank you.

• Does this answer your question? Solving coupled differential equations in Python, 2nd order – Bill Greene Jan 29 at 14:12
• @BillGreene Thank you for your reply : unfortunately, I have already seen this conversation but it's not the same kind of simplification that is needed... I think I can apply a shooting method algorithm on this type of problem but I don't see how to simplify it. Thank you again. – Wiss Jan 29 at 14:21
• So this is a boundary value rather than an initial value problem? In any case, you should edit your question to include the boundary or initial conditions. What do you mean by "may need a simplification that I don't see"? Converting higher order ODE to first order form is something that has been discussed widely. I strongly suggest you do not try to implement your own algorithm for solving this system unless there is some feature that is not evident from your post. – Bill Greene Jan 29 at 14:44
• @BillGreene Yes it is a Boundary value problem : I have updated my post in order to clarify the boundary conditions. I mean that maybe I need a transformation to reduce the order of each equation in order to simplify it. In fact I used to solve linear BVP by a shooting method algorithm so I have already done it before but this particular system doesn't allow me to apply the shooting method so I am a little bit lost in order to find a strategy to solve it. Thank you. – Wiss Jan 29 at 16:25
• Hello again @LutzLehmann : the x is just the multiplication sign. The Pr has to be superior than one so I attend to do a simulation with Pr varying between 1 and 10. Thank you – Wiss Jan 30 at 10:02

There is no higher magic necessary, just transcribe into the canonical first-order system, encode the boundary conditions, make a reasonable initial guess of the solution shape and call the BVP solver

Pr = 5

def odesys(t,u):
F,dF,ddF,θ,dθ = u
return [dF, ddF, θ-0.25/Pr*(2*dF*dF-3*F*ddF), dθ, 0.75*F*dθ]

def bcs(u0,u1): return [u0[0], u0[1], u1[2]-1, u0[3]-1, u1[3]]

x = np.linspace(0,1,4)
u = [0.5*x*x, x, 0*x+1, 1-x, 0*x-1]

res = solve_bvp(odesys,bcs,x,u, tol=1e-5)
print(res.message)

Then plotting the solution gives

• Once again you allowed me to solve a system of differential equations... In fact I used to solve BVP with my own code that works fine with simple cases but for this one I had some problems with the boundary conditions... I have to work on that for the future : Thank you again – Wiss Jan 30 at 11:07

What does differentiating the first equation once to give $$\theta'$$ and twice for $$\theta''$$ and plugging into the second equation to give a single equation for $$F$$ give? Seems like solving a single equation for $$F$$ might be your best approach (or vice versa).

• Thank you for your answer : I tried this option but it ended with a differential equation of order 5 and I don't have the required number of boundary values in order to solve it. – Wiss Jan 29 at 18:04
• Seems like you have plenty of BCs, but some of them will end up being of mixed type. – Bill Barth Jan 29 at 18:10