# 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 Commented Jan 29, 2021 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
Commented Jan 29, 2021 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. Commented Jan 29, 2021 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
Commented Jan 29, 2021 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
Commented Jan 30, 2021 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
Commented Jan 30, 2021 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
Commented Jan 29, 2021 at 18:04
• Seems like you have plenty of BCs, but some of them will end up being of mixed type. Commented Jan 29, 2021 at 18:10