integrating of 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})