# 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- Ignore the multiple plots, because of the parameter, $$P$$, but the shape should match. Posted on Categories Articles