differential equations - (NDSolve) Non-linear 2nd order ODE, regular singular point (looking for good methods for this problem)

I am solving this set of non-linear 2nd order ODE by NDSolve, $$r^2\frac{d^2f}{dr^2} = 2f(1-f)(1-2f)+\frac{r^2}{4}h^2(f-1)$$ $$\frac{d}{dr}\left[r^2\frac{dh}{dr} \right]=2h(1-f)^2+\lambda r^2(h^2-1)h$$ $$ f(0)=h(0)=0, \quad f(\infty)=h(\infty)=1$$

From the series solution, I also know the behavior at $r\approx0$, $$ f\approx \alpha r^2, \quad h\approx \beta r$$ Therefore the singularity at $r=0$ is a regular singular point.

I tried using NDSolve,

\[Lambda] = 625/2048;
eqn = {r^2*f''[r] == 
2 f[r] (1 - f[r]) (1 - 2 f[r]) - r^2/4 (h[r])^2 (1 - f[r]), 
D[r^2*h'[r], r] == 2 h[r] (1 - f[r])^2 + \[Lambda]*r^2 (h[r]^2 - 1) h[r]};
bc = {f[0.000001] == 0, f[20] == 1, h[0.000001] == 0, h[20] == 1};
sol = NDSolve[{eqn, bc}, {h, f}, {r, 0.000001, 20}]

However, it doesn't solve it says At r =...., step size is effectively zero; singularity or stiff system suspected.

I also tried,

NDSolve[{eqn, bc} /. {f -> (#^2 g[#] &), h -> (# j[#] &)}, {g, j}, {r, r1, r2}]

with boundary condition at r1 and r2. Still it doesn't give me the solution.

I guess there are some subtleties about non-linear ODEs. Since it is non-linear, the solution is highly sensitive to boundary conditions. Unless I set the exact boundary condition, it won't give me the right answer. However, the ODE above has a regular singular singular point at $r=0$, which I don't know how to impose the exact boundary condition at $r=0$ in mathematica. It stops calculation immediately if I set the boundary condition at $r=0$.

I am looking for help from you guys, especially those had the same problem before. My sincere thank for your help!

Do you know if the method of relaxation help in this case?

Update: The method of relaxation works well for this type of problem, however it seems that Mathematica does not implement this for NDSolve algorithm.