a := evalf(Constant(hbar)*Constant(m[P])*(8*Pi*(1/2))^(1/2)/Constant(c));

#mp := (8.5*1.783)*10^(-33);

chi := a/mp^2;

pctom := 3.08568025*10^16;

solmass := 1.98892*10^30;

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;

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

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;

bettermcproc := proc (R, Theta0, Beta, ThetaR) local rcd; Wsol(parameters = ['theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR]); rcd := [fsolve(diff(safervc(r, Theta0, Beta, ThetaR), r), r, r = 1/1000 .. 60)]; if not type(rcd, list(numeric)) then rcd := [fsolve(proc (R) options operator, arrow; evalf((D(safervc(R, Theta0, Beta, ThetaR)))(R)) end proc, 1/1000 .. 60)]; if not type(rcd, list(numeric)) then error "fsolve stage failed for parameter values %1", Msol(parameters) else Msol(rcd[1]) end if else Msol(rcd[1]) end if end proc;

rcproc := proc (R, Theta0, Beta, ThetaR) Wsol(parameters = ['theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR]); fsolve(diff(safervc(r, Theta0, Beta, ThetaR), r), r, r = 1/1000 .. 60) end proc;

#rhproc := proc(R, Theta0, Beta, ThetaR) ... end; # Looks like there's a problem here

mhproc := proc (R, Theta0, Beta, ThetaR) local rcd; Wsol(parameters = ['theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR]); rcd := [fsolve(diff(safervc(r, Theta0, Beta, ThetaR), r), r, r = 1000000 .. 1000000000)]; if not type(rcd, list(numeric)) then rcd := [fsolve(proc (R) options operator, arrow; evalf((D(safervc(R, Theta0, Beta, ThetaR)))(R)) end proc, 1000000 .. 10000000)]; if not type(rcd, list(numeric)) then error "fsolve stage failed for parameter values %1", Msol(parameters) else Msol(rcd[1]) end if else Msol(rcd[1]) end if end proc;

vc:=eval(sqrt(M(r)/r),dsol):

safervc:=proc (R, Theta0, Beta, ThetaR) Wsol(parameters = ['theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR]); vc(R) end proc ;