How to solve a system of ODEs with BCs, without bvp4c?

Hi, I solved a system of ODEs with BCs, by using bvp4c.

The system is a bit complicated, consisting of 21 non linear equations, with 21 variables depending on space.

This is the reason why the simulated solution is extremely dependant on the initial guess, with an high maximum residual compared to requested accuracy. This is very critical, because I don't know with precision what values should I enter as initial guess in order to have a "good" solution.

Hence, I'm looking for a strongest method to solve my sistem.

Do you suggest me something else, if it exists? Is it possible by using a standard ode solver?

Thank you!