# 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}$$ The boundary conditions of the problem are-

At $$x=0$$, $$v=w=0$$ and $$\lambda_2=\lambda_0$$ where $$\lambda_0$$ is some constant.

End the iteration when $$\lambda_2=1$$ and check value of $$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]][]
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- Ignore the multiple plots, because of the parameter, $$P$$, but the shape should match. I am getting a plot like this- • The initial conditions matter a lot. Try setting the start value to where you specify the initial conditions to: 0.001, 0.0009,0.0008, 0.0007, .. and you will see that the system is sensitive to the initial conditions Oct 12 '20 at 15:56
• That was a good advice. I tried that, but of no use. Also, I should have specified the initial conditions in the question itself. I am editing the question rn.
– user75147
Oct 12 '20 at 18:13
• I think I know what you mean, but your question is horribly formulated. If you first make a calculation with a start value of 0.1, then 0.01, then 0.001 e.t.c you will see that the solution diverges. Therefore, you can NOT start with 0. Oct 12 '20 at 18:44
• Yes, I agree. I could have put up the problem a bit more clearer. As per your advice, I am trying to see the solution by slowly increasing the initial conditions.
– user75147
Oct 12 '20 at 19:47
• " it should be the step-size, I think." I don't think so. The ODE solver of NDSolve is quite robust and should always be the last thing to tackle. Better to double check the system itself. Oct 13 '20 at 2:27

The problem is ill-posed; $$v$$ cannot vanish at $$x = 0$$.

Assuming that $$v(0) = 0$$, $$w(0) = 0$$, and $$\lambda_2(0) = \lambda_0 = 1.5$$, we have $$v'(x) \sim \frac{1.08306}{x}$$ for small $$x$$, which is inconsistent with $$v(0) = 0$$.

derivativeList = {
F1[x, v[x], w[x], lambda2[x], P],
F2[x, v[x], w[x], lambda2[x], P],
v[x]/x
};

initialRules = {
v[x] -> 0,
w[x] -> 0,
lambda2[x] -> lambda0
};

derivativeList /. initialRules // FullSimplify
(* {0. + 1.08306/x, 0.308318, 0} *)

• Thanks for the information. The system of equations coming out has a problem. I have to rethink it.
– user75147
Oct 18 '20 at 4:47