Question: Using a numerical solution on another function

Hello all

I have a set of differential equations and its numerical solution, for two functions. Here's the worksheet (just a part of it, I don't need the rest)

##########################################

eq1 := diff(W(r), r) = -(1-beta*(W(r)-W0))*(M(r)+4*Pi*r^3*p(r))/(beta*r*(r-2*M(r)));

eq2 := diff(M(r), r) = 4*Pi*r^2*rho(r);

W0 := solve(thetaR = theta0-W0, W0);

r0:=10^(-5);

rf:=10^20;

# ini := evalf({M(r0) = 0, W(r0) = W0});
ini := [M(r0) = 0, W(r0) = W0];

thetadef := theta(r) = W(r)+thetaR;

#nuvdef:=nuv(z)=z^(-2*alpha)*ApproximateInt(r^(2*alpha-2)*3*l*(4*Pi*rh^3)^(-1)*(1+r^2/rh^2)^(-5/2)*M(r),r=z..10^5,output=sum,partition=50):

pdef := p(r) = subs(thetadef, value((8*sqrt(2)*Pi*(1/3))*(beta/(1-beta*(W(r)-W0)))^(5/2)*ApproximateInt((1+(1/2)*beta*x/(1-beta*(W(r)-W0)))^(3/2)*x^(3/2)*(1-exp(x-W(r)))/(exp(x-theta(r))+1), x = 0 .. W(r), output = sum, partition = 50))):

rhodef := rho(r) = subs(thetadef, value(4*sqrt(2)*Pi*(beta/(1-beta*(W(r)-W0)))^(3/2)*ApproximateInt((1+(1/2)*beta*x/(1-beta*(W(r)-W0)))^(1/2)*(1+beta*x/(1-beta*(W(r)-W0)))^2*x^(1/2)*(1-exp(x-W(r)))/(exp(x-theta(r))+1), x = 0 .. W(r), output = sum, partition = 50))):

# subs(pdef, rhodef, [eq1, eq2, op(ini)]): # sometimes the order of the substitutions can matter

# dsys := evalf(subs(pdef, rhodef, nuvdef, {eq1, eq2, op(ini)}));
dsys := subs(pdef, rhodef, {eq1, eq2, op(ini)}):
dsol:=dsolve( dsys
, 'type' = numeric
, 'output' = listprocedure
, 'range' = r0 .. rf
, 'parameters' = [theta0,beta,thetaR]
 , 'complex' = false
# , 'compile' = true # these options sometimes help increase efficiency
) ;

Wsol:=eval(W(r),dsol):

Msol:=eval(M(r),dsol):

thetasol := subs(W(r) = Wsol(r), rhs(thetadef));

psol:=subs(W(r)=Wsol(r),rhs(pdef)):

rhosol:=subs(W(r)=Wsol(r),rhs(rhodef)):

betterpsol := proc (R, Theta0, Beta, ThetaR) Wsol(parameters = ['theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR]); evalf[15](eval(psol, [r = R, 'theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR])) end proc;

betterrhosol := proc (R, Theta0, Beta, ThetaR) Wsol(parameters = ['theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR]); evalf[15](eval(rhosol, [r = R, 'theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR])) end proc;

##############################

Now, I need to define another function that depends on the solution for M(r), like

f:=int(h(r,parameters)*M(r), r=z..infinity)

where h(r) depends on different parameters than the system.

Then I need another function

y:=G(R)*int(b(z,same parameters as h)*f(z))

What I would like is to be able to work with this function by calling it like y(R,theta0,beta,thetaR, 'extra parameters'), or something similar. To explain better, I don't want to define the extra parameters early, but only when calling the function y(R), similar to betterpsol and betterrhosol.

I tried a couple of things to implement this, but it didn't work and I'm out of ideas.

Thanks for the help