bvpinit giving wrong behavior!

I'm tryping to solve the following nonlinear Diff. Eq:

, using the function bvp4c and the boundary conditions: . After finding α, I need to use and to model something. Also, t is used to help simplify the math, but it is related to distance r via: . Here, is just some characteristic length scale for the system.

The issue is I expect linear behavior for around (being ) since α will be small here. This can be seen by realizing continuity of α and our boundary condition on the left implies this. Indeed, our differential equation reduces to (near this point):

The plot I have included (r ) illustrates the behavior is quadartic near ... I can not figure out exactly why!

I included my code and setup as well. The only reason I can see something may go wrong is I mapped to using tanh. I tried running bvp4c with something big, so I could still treat it numerically, but no luck... I kept getting a singular Jacobian.

Any help will be greatly appreciated!