# Integrating the second derivative of numerical solution

pde1 = -y1''[x] - (2*y1'[x])/x + ((y1[x])^3 + y2[x])y1[x] == 0;
pde2 = y2''[x] + (2y2'[x])/x - (y1[x])^3 == 0;
sol = NDSolve[ {pde1, pde2, y1[1] == 0.001, y2[1] == -0.001,
y1'[0.001] == 0.001, y2'[0.001] == 0.001}, {y1, y2}, {x,0.001, 20}]


I need to plot the values of Integrate[y1''[x] x^2, {x, 0.001, 20}].

Try NIntegrate

NIntegrate[y1''[x] x^2 /. sol[[1]], {x, 0.001, 20}]
(*-0.000826417*)

• Thank you so much Nov 18 '20 at 16:53

You can integrate by parts twice and get an expression for the antiderivative:

parts[u_, v_, {x_, n_}] :=
Sum[(-1)^m D[u, {x, m}] Nest[Integrate[#, x] &, v, m + 1], {m, 0,
n - 1}] + (-1)^n Integrate[
D[u, {x, n}] Nest[Integrate[#, x] &, v, n], x];

parts[x^2, y1''[x], {x, 2}]
int1[x_] = % /. First[sol];
int1[20] - int1[0.001]

(*
2 Integrate[y1[x], x] - 2 x y1[x] +  x^2 y1'[x]

-0.000826417
*)


Plot:

Plot[int1[x], {x, 0.001, 20}]

• Thank you so much Nov 18 '20 at 16:53
• @BahiMido You're welcome. :) Nov 19 '20 at 5:25