You can add a extra differential equation f'[x]==y2[x]*y1[x]^3* x^2
in your equations and solve f'[x]
at the same time.
so we need to set the initial value f[0.001]=c
and c
is a parameters,so we use ParametricNDSolve
Clear[pde1, pde2, pde3, sol];
pde1 = -y1''[x] - (2*y1'[x])/x + ((y1[x])^3 + y2[x]) y1[x] == 0;
pde2 = y2''[x] + (2 y2'[x])/x - (y1[x])^3 == 0;
pde3 = f'[x] == y2[x]*y1[x]^3*x^2;
sol = ParametricNDSolve[{pde1, pde2, pde3, y1[1] == 0.001,
y2[1] == -0.001, y1'[0.001] == 0.001, y2'[0.001] == 0.001,
f[0] == c}, {y1, y2, f}, {x, 0.001, 20}, {c}]
f[0][20] - f[0][0.001] /. sol
f[1][20] - f[1][0.001] /. sol
Plot[f[0][x] /. sol, {x, 0.001, 20}]
Plot[f[1][x] /. sol, {x, 0.001, 20}]
-2.38247*10^-9
-2.3811*10^-9
{fy1, fy2} = Values@First@sol; NIntegrate[fy2[x]*fy1[x]^3 x^2, {x, 0.001, 20}]
.What you're asking for gives a single number so I don't know why you're asking for a plot. $\endgroup$ – flinty Dec 2 '20 at 18:52