Hello

This is actually somewhat the same question I asked some months ago. I found a way to go around it, but now unfortunately I'm back to the same problem, and I didn't have any answers when I asked.

Here it is a system that I solve numerically:

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

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 := [M(r0) = 0, W(r0) = W0];

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

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))):

dsys := subs(pdef, rhodef, {eq1, eq2, op(ini)}):

dsol:=dsolve( dsys

, 'type' = numeric

, 'output' = listprocedure

, 'range' = r0 .. rf

, 'parameters' = [theta0,beta,thetaR]

, 'complex' = false

) ;

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)):

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

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 a new variable that obeys this equation:

diff(v(r)*nu(r),r)=-M(r)*nu(r)/r^2, with nu(r) a known function

with a boundary condition v(infinity)=0. I don't know the value of this function at r0, so when I tried this first I just added another initial condition (vrad(10^20)=0) but it didn't work. However, I know the solution for v(r) in terms of an integral on M(r), the same M(r) of the original system:

v(R)=(nu(R))^(-1)*int(nu(r)*M(r)/r^2,r=R..infinity).

I tried to do this as a procedure

v:=proc(X,Theta0,Beta,ThetaR,extra parameters from nu)

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

evalf(eval(int...,r=x..infinity,[x = X, 'theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR,parameters from nu(r)]))

end proc:

This takes an unreasonable amount of time (more than half an hour, at which point I aborted). Did I have to declare 'r' as a local variable in the procedure?

I also tried a shooting solution, including v(r0)=alpha on ini, and bc:={v(10^20)=0}, and the equation for v(r) in the system, then

shoot(dsys,[M(r),W(r),v(r)],bc,ini,[alpha=1e-2],output=listprocedures,range=r0..rf,parameters=[same+extra ones]);

I wanted to call v(r) the same way I call betterrhosol and betterpsol, like v(r,theta0,beta,thetar,extra parameters).

I am now running it and it hasn't finished yet after half an hour. Any help?

Thanks a lot