How can I "simplify" the resolution of this system of odes?

Hi, as you can see in the attached file, I've a system of second order ODEs with BCs.

I tried to solve it with bvp4c, firstly converting second order equations in first order. So, I don't ask you how to solve it with bvp4c, since I've already did it. My biggest problem is that final solution is extremely dependent on the initial guess of solution since the equations are non linear, and numerical uncertainty is too high.

How can I find a way to make the resolution of the system "easier", for example solving the equations through a discretization?

Thank you!