Boundary value problem for natural convection by shooting method

Question

I am trying to resolve a special case of boundary value problem which is the boundary layer's equations for natural convection :

Thanks to a contributor of this forum @LutzLehmann, this set of equations is solved by using the function solve_bvp from the scipy library.

Unfortunately, I don't obtain the right result with these particular boundary conditions :

It seems that the initial guess for the solver doesn't work in this case and I wonder what could suit better.

Here is the code :

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt


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[1], u0[3]-1, u1[3]]

x = np.linspace(0,8,25)
u = [x*x, np.exp(-x), 0*x+1, 1-x, 0*x-1]

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

print(res.message)

plt.subplot(2,1,1)
plt.plot(res.x,res.y[3], color='#801010', label='$\Delta T$')
plt.legend()
plt.grid()

plt.subplot(2,1,2)
plt.plot(res.x,res.y[1], '-', color='C0', label="$F'$")
plt.legend()
plt.grid()

And here are the wrong plots obtained :

Could someone help me by suggesting a better initial case for this problem ?

Thank you for your help,

PS : For the proposition :

u = [0.5*x*x*np.exp(-x), x*np.exp(-x), np.exp(-x), np.exp(-x), -np.exp(-x)]

I got this :

Compared to the literature, I am expecting this result (where g is here theta in the set of equations) :