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 ;