@Carl Love Thanks, this worked. But what if the function has more than one parameter, like F(x,a,b,c)? zip does not seem to work then

@Carl Love 

wouldn't this work only if the value of the parameter a starts with 1 and is an integer? If I change to

xlist:=[seq(X1(i), i=1..10,0.5)]:

it doesn't work anymore

The maxfun=0 option works fine, thanks a lot. I tried to use a large number but still wasn't computing. 

As for your procedural solutions, the first one doesn't work. I tried tweaking a bit, but I still get the same error: 

Error, (in f) unable to store 'HFloat(22509.918246846857)*L' when datatype=float[8]

The second one takes a very long time to calculate, and I'm not really sure I'm doing it right. I called the solution in a point as res(10) for example, but it took a very long time.

I tried to call ysol:=eval(y(x),res(2.3)) for example, but it doesn't get any better. ysol(10) takes more than 5 minutes, the same with ysol:=rhs(res₂). 

Thanks for the help

The maxfun=0 option works fine, thanks a lot. I tried to use a large number but still wasn't computing. 

As for your procedural solutions, the first one doesn't work. I tried tweaking a bit, but I still get the same error: 

Error, (in f) unable to store 'HFloat(22509.918246846857)*L' when datatype=float[8]

The second one takes a very long time to calculate, and I'm not really sure I'm doing it right. I called the solution in a point as res(10) for example, but it took a very long time.

I tried to call ysol:=eval(y(x),res(2.3)) for example, but it doesn't get any better. ysol(10) takes more than 5 minutes, the same with ysol:=rhs(res₂). 

Thanks for the help

You're right, I didn't realize that the vertical axis had logarithmic ticks for some reason. Thank you

You're right, I didn't realize that the vertical axis had logarithmic ticks for some reason. Thank you

The definition of r0 and rf are missing, I usually use

r0:=10^(-4)

rf:=10^(17)

The definition of r0 and rf are missing, I usually use

r0:=10^(-4)

rf:=10^(17)

@acer Well, since I got no answer to those questions I thought I'd try again, but making it simpler. And yes, it is a followup from april. 

The variable in M(r) doesn't matter, it's a dummy variable anyway. But when I tried M(s) or something similar it gave me an error, since this is the same M as defined in eq2. 

As for ApproximateInt: i was trying with evalf(int(...)) before, but it was taking a loooooong time, maybe because the upper limit is a function to be solved in the system. Someone suggested me ApproximateInt, and so far I have no complaints. But I got a better pc now, do you think evalf(int(..)) will give me a more accurate result? 

And for the problem itself, any ideas?

Thanks for the answer, I was usind the "read" command. Anyways, I solved it, simple mistake of my part.

Thanks for the answer, I was usind the "read" command. Anyways, I solved it, simple mistake of my part.

Forgot to add the input.txt, here it is

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 ;

Download input.txt

@PatrickT 

Thanks for the help, but how can it be a file permission problem? I tried on a couple other computers and it keeps failing. I tried just one value

bettermcproc(10,2.4e5,10e-6,-60)

with the input and it fails; it works fine if I don't use the input file. Well, I'll try opening another post for this

Thanks a lot for the help

@PatrickT 

Thanks for the help, but how can it be a file permission problem? I tried on a couple other computers and it keeps failing. I tried just one value

bettermcproc(10,2.4e5,10e-6,-60)

with the input and it fails; it works fine if I don't use the input file. Well, I'll try opening another post for this

Thanks a lot for the help

@PatrickT 

well, the .txt is the same as you wrote before, I'll copy it again just in case

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

) ;

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 .. 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;

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 ;

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

main file

restart;
### Load Packages
with(ScientificConstants):
with(Student:-Calculus1,ApproximateInt);

with(plots):
read "C:/maple/testing/input.txt" : #the path on my file is correct, I just forgot and the file is in the university pc

mcseq:=[seq(bettermcproc(10,i,10^(-6),-60),i=2*10^5..3*10^5,10^3)]:

Then it fails, giving me the error message from bettermcproc. The point is, I have another worksheet where I solve the equations and have the procedures in the worksheet itself, and there it works fine. Apparently there's some problem for high values of that parameter, but I don't know why