@bernie

It looks like your input file is mainly procedures, so I think this should work. I have only tested it on plot p1.

Remarks:

I prefer to load packages at the top, to keep track, in case of interferences (can happen), also I prefer to load only specific commands if I'm using only one or two from a package (see below where I load ApproximateInt only). Should there be an interference, it's easier to test if you know exactly which packages are loaded, especially when you load external files and the commands become "hidden" from view.

I prefer to enter dsolve inputs as lists, rather than sets, to keep the ordering consistent across sessions. If you're going to use evalf, I suggest you do it at the last moment. I don't think you need to evalf here, seemed to work without it.

Do scroll down, I have written one or two comments in passing.

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

restart;

### Load Packages

with(ScientificConstants):

with(Student:-Calculus1,ApproximateInt);

### Load my procedures

read "C:/maple/testing/input.txt" :

p1 := plots:-semilogplot(thetasol1(r, 6*10^4, 10^(-6), -50), r = r0 .. 10, axes = boxed, color = black, thickness = 2):

p1; # works okay here

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

The following is my input.text, which is just a text file.

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

Wsol(parameters = [theta0 = 24, beta = 3.17*10^(-7), thetaR = -31]); current := convert(Wsol(parameters), `global`); evalf(Wsol(.645*10^8)); evalf(Msol(.645*10^8)); evalf(Constant(c)^2*chi*Msol(.645*10^8)/(Constant(G)*solmass));

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 .. 10)]; 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 .. 10)]; 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 .. 10) 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;