@Preben Alsholm I actually made it work like this:

nu:=r->f(r) #the function doesn't matter here

nuvrad:=evalf(int(nu(s)*Msol,s=x..10^20)):

vradsol:=proc(X,Theta0,Beta,ThetaR,L,Rh) #last ones are the parameters in f(r)

Wsol(parameters = ['theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR]);

evalf(eval(1/nu(X)*vrad(X),[x=X, 'theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR, 'l'=L, 'rh'=Rh]))

end proc:

and it works fine. However, I need again an integral of this function vradsol. I managed to do it the same way, but I simply can't plot the function, it takes too long. Any suggestions?

Thanks to everyone, all of these solutions worked! I'd rather the exponent didn't have brackets, so for my purposes I think acer's solution was the best. Didn't know anything this simple could get complicated. Aren't there any plans for the maple team to give us more control over the plot legend, so we can place it anywhere we want?

@acer The variable of integration in vrad(x) is r, the same as in dsolve,numerical. This is because I need to integrate M(r), one of the solutions of the system.

I tried writing an expanded ode system, but it doesn't work it says that it cannot solve a parametric ode system. Could it be because the boundary condition for vrad is vrad(infinity)=0? (I actually used a number, not infinity, still didn't work).

As for the accuracy, I just assumed that ApproximateInt was less acurate than int, and that I was pretty happy with the accuracy I had