I have figured out the answer to my own question. Starting with the numerical solution, which in my case looks like

[y = -3., U(y) = 0.395539692200000e-4, W(y) = -5.90967381270470, phi[1](y) = 1.13186853677760, phi[2](y) = 0.]

and is called 'sol', we can define W from

solW := r->rhs(sol(r)[3]);

Then we can plot and integrate using 

plot(solW,ymin..ymax)

and 

g := proc (k)
evalf(Int(y -> solW(y), ymin .. k))
end proc;

g(ymax); # THIS IS THE INTEGRAL

Notice that Int is capatalized!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!