for exemple i have pde1 = -y1”[x] – (2*y1′[x])/x + (10^-19*(y1[x])^2 + y2[x])*y1[x] == 0;
pde2 = y2”[x] + (2*y2′[x])/x – (10^-25)(y1[x])^2 == 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, 30}]
{{y1 -> InterpolatingFunction[{{0., 30}}, <>],
y2-> InterpolatingFunction[{{0., 30}}, <>]}}
I need to plot the integration of ( solution(y1) x^2)