Problem in getting correct plot to the numerical solution of a non-linear differential equation in Mathematica


I have a system of equations to be solved numerically-
$$
v'(x) = F_1(v,w,x,lambda_2,P) \
w'(x) = F_2(v,w,x,lambda_2,P) \
lambda_2′(x) = frac{v}{x}
$$

I have found out the system of equations symbolically. I am attaching the code-

    lambda1 = Sqrt((lambda2(x) + x lambda2'(x))^2 + (eta'(x))^2);
    lambda3 = 1/(lambda1 lambda2(x));
    I1 = lambda1^2 + lambda2(x)^2 + lambda3^2;
    I2 = 1/lambda1^2 + 1/lambda2(x)^2 + 1/lambda3^2;
    alpha = 0;
    W = -(I1 - 3) - alpha (I2 - 3) + 0.5 x lambda2(x)^2 eta'(x) P;
    eq1 = Numerator(Together(D(W, lambda2(x)) - D(D(W, lambda2'(x)), x)));
    eq2 = Numerator(Together(D(D(W, eta'(x)), x)));
    eq1mod = Expand(
       eq1 /. {eta'(x) -> w, eta''(x) -> w', lambda2'(x) -> v/x, 
         lambda2''(x) -> (v' - v/x)/x, lambda2(x) -> lambda2});
    eq2mod = Numerator(
       Together(
        Expand(eq2 /. {eta'(x) -> w, eta''(x) -> w', lambda2'(x) -> v/x, 
           lambda2''(x) -> (v' - v/x)/x, lambda2(x) -> lambda2})));
    f1 = Coefficient(eq1mod, v');
    f2 = Coefficient(eq1mod, w');
    f3 = Simplify(-(eq1mod - f1 v' - f2 w'));
    f4 = Coefficient(eq2mod, v');
    f5 = Coefficient(eq2mod, w');
    f6 = Simplify(-(eq2mod  - f4 v' - f5 w'));
    F1sym = -(f2 f6 - f3 f5)/(f1 f5 - f2 f4);
    F2sym = (f1 f6 - f3 f4)/(f1 f5 - f2 f4);

I know I am correct upto this point. Now I converted the 2 symbolic expressions to functions F1(x,v,w,lambda2,P) and F2(x,v,w,lambda2,P).

    exprToFunction(expr_, vars_) := 
      ToExpression(
       ToString(FullForm(expr) /. MapIndexed(#1 -> Slot @@ #2 &, vars)) <>
         "&");
    F1 = exprToFunction(F1sym, {x, v, w, lambda2, P});
    F2 = exprToFunction(F2sym, {x, v, w, lambda2, P});

On closer look, I found the system to be singular at x=0. So I started the calculation from x=0.001. I am putting the code below-

    P = 0.5;
    lambda0 = 1.5;
    sol = NDSolve({
        v'(x) == F1(x, v(x), w(x), lambda2(x), P),
        w'(x) == F2(x, v(x), w(x), lambda2(x), P),
        lambda2'(x) == v(x)/x,
        v(0.001) == 0, w(0.001) == 0, lambda2(0.001) == lambda0,
        WhenEvent(lambda2(x) == 1, "StopIntegration")},
       {v, w, lambda2}, {x, 0.001, 20});
    end = InterpolatingFunctionDomain(First(lambda2 /. sol))((1))
    Plot(Evaluate(lambda2(x) /. sol), {x, 0.001, 1}, Frame -> True)

Now, after solving, I find the output plot not correct. I don’t know the real reason, but I think, in the first step itself, there is a wrong step-size selected. How to correct this problem?
The correct image should be this-
enter image description here

Ignore the multiple plots, because of the parameter, $P$, but the shape should match.