I am trying to solve a system of linear ODEs with two unknown parameters using Python. But I am facing some problems. Below is my code:

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

#parameters

chi = 1; #ratio of reference density at two boundaries
m = 1.495; #polytropic constant
th = -chi**(-1/m)*(-1+chi**(1/m)); #degree of compressibility
Nrho = np.log(abs(chi)); #the number of density scale heights

l = 2; #horizontal wavenumber along x
k = 0; #horizontal wavenumber along z
a = np.sqrt(l**2+k**2);

Pr = 0.1;
Ta = 1e05;
ph = np.pi/4;


def fun(z,w,p):
    lm = p[0];
    Ra = p[1];
    return np.vstack((w[1],w[2],w[3],w[5],w[7],-1j*Pr*Ta**(1/2)*(l*(1+z*th)*np.cos(ph)-1j*m*th*np.sin(ph))/(Pr*(1+z*th))*w[0]-\
                      Pr*Ta**(1/2)*(1+z*th)*np.sin(ph)*w[1]/(Pr*(1+z*th))+(a**2*Pr*(1+z*th)+lm+z*th*lm)/(Pr*(1+z*th))*w[4]+m*Pr*th*w[5]/(Pr*(1+z*th)),\
                      a**2*Ra*w[6]+(2*a**2*m**2*th**2/(3*(1+z*th)**2)-(3*(m-2)*m*th**4+a**4*(1+z*th)**4)/(1+z*th)**4+lm/Pr*(-a**2-m*th**2/(1+z*th)**2)+\
                        1j*k*m*Ta**(1/2)*th*np.cos(ph)/(1+z*th))*w[0]+(3*(m-2)*m*th**3/(1+z*th)**3+2*a**2*m*th/(1+z*th)+m*th*lm/(Pr*(1+z*th)))*w[1]+\
                      (2*a**2+(4-m)*m*th**2/(1+z*th)**2+lm/Pr)*w[2]-2*m*th*w[3]/(1+z*th)+1j*l*Ta**(1/2)*np.cos(ph)*w[4]+Ta**(1/2)*np.sin(ph)*w[5]\
                      ,(a**2+(1+z*th)**m*lm)*w[6]-th*w[7]/(1+z*th)-(1+z*th)**(-1+m)*w[0]));


def bc(wa,wb,p):
    lm = p[0];
    Ra = p[1];
    return np.array([wa[0],wb[0],wa[2]+m*th/(1+th*0)*wa[1],wb[2]+m*th/(1+th*1)*wb[1],wa[5],wb[5],wa[6],wb[6],wa[1]-1,wa[7]-1],dtype=complex);


z = np.linspace(0,1,1000);
w = np.zeros((8,z.size), dtype=complex);

w[0] = np.sin(np.pi*z)

sol = solve_bvp(fun, bc, z, w, p =[8+11j,1e05]);

if sol.status != 0:
    print("WARNING: sol.status is %d" % sol.status)
print(sol.message)
print(sol.p[0],sol.p[1])

The problems here: 1.) I don't understand how to get the algorithm to converge. 2.) How to set the initial guesses? Please help me. Using Mathematica I was able to find lm=8.0489 + 11.3672*I and Ra=2e05.