@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

@PatrickT 

Ok, it worked perfectly until I tried to do

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

and it fails. On the worksheet I had (without the .txt input) this works fine. Apparently it doesnt work for high values of theta, but I don't know why.

@PatrickT 

Ok, it worked perfectly until I tried to do

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

and it fails. On the worksheet I had (without the .txt input) this works fine. Apparently it doesnt work for high values of theta, but I don't know why.

So far it works perfectly, thanks a lot. You've saved me a lot of time and work.

Thanks

So far it works perfectly, thanks a lot. You've saved me a lot of time and work.

Thanks

Ok, here's my worksheet

scalingedited.mw    It has a lot of stuff, and I would like to make it neater by splitting the things I've done in different worsheets without having to solve the differential equations in the beginning.

I need to learn to post my questions more clearly...what I would like is to have a file that I could read or import straight away, that looks like this:

input.mw 

So I could somehow use all this information in other worksheets. Hope it's clear now

Tried to post the second one as comment and got some error, so there's a link instead

Thanks

Maybe I didn't explain myself right. The data is fine, I want to use all the data. The point is, the model function f(x,parameters) depend on a variable x; it's this variable that I want to set the range. I realize now that using xdata may have confused things a bit. For this problem, small values of x in the function don't matter, I want Fit to use the variable x above a certain value. 

If I don't get how Fit works and this can't be done, maybe there's a way out. The problem is that for small values of x, f(x) is decreasing, while the data is increasing. Only after a certain point is that f(x) is increasing again. If I just try Fit, it gives me complex values or it says that the initial values are wrong; however I know for sure that there's a fit, no matter how bad it is. I'm at home now and very smartly forgot the file in university, tomorrow I can post the pointplot and a plot I made guessing some values for the parameters, and they are similar enough to give me a fit, I think.

Maybe I didn't explain myself right. The data is fine, I want to use all the data. The point is, the model function f(x,parameters) depend on a variable x; it's this variable that I want to set the range. I realize now that using xdata may have confused things a bit. For this problem, small values of x in the function don't matter, I want Fit to use the variable x above a certain value. 

If I don't get how Fit works and this can't be done, maybe there's a way out. The problem is that for small values of x, f(x) is decreasing, while the data is increasing. Only after a certain point is that f(x) is increasing again. If I just try Fit, it gives me complex values or it says that the initial values are wrong; however I know for sure that there's a fit, no matter how bad it is. I'm at home now and very smartly forgot the file in university, tomorrow I can post the pointplot and a plot I made guessing some values for the parameters, and they are similar enough to give me a fit, I think.

@pagan Thanks a lot, I think I understand now how to do it. To answer your questions:

I'm not an expert maple user, this is the first time I'm using it for something more complicated than just solving ODEs numerically. So, if I remember correctly, I had some problems with evalf(int(...)), maybe because the upper limit is actually a solution of the system, so someone told me to use ApproximateInt instead. As for the proc(r) parts, I copied it from another worksheet of a guy doing similar stuff. Now that I know, I'll start using unnaply.

And to clarify, I meant quicker for me, not for the computer. Without the loops, I'd have to input the parameters everytime I needed a different solution, now I can just let the machine work it out.

Thanks a lot for the help

@pagan Thanks a lot, I think I understand now how to do it. To answer your questions:

I'm not an expert maple user, this is the first time I'm using it for something more complicated than just solving ODEs numerically. So, if I remember correctly, I had some problems with evalf(int(...)), maybe because the upper limit is actually a solution of the system, so someone told me to use ApproximateInt instead. As for the proc(r) parts, I copied it from another worksheet of a guy doing similar stuff. Now that I know, I'll start using unnaply.

And to clarify, I meant quicker for me, not for the computer. Without the loops, I'd have to input the parameters everytime I needed a different solution, now I can just let the machine work it out.

Thanks a lot for the help

Ok, I can make the loops, but as seen on the worksheet, it doesn't return different values for different parameters. Also, for some reason this values are wrong. I did just for rc as you can see, but this happens to every quantity. For p0 and rho0 I don't even get numerical values. I interrupted the computation because it was taking too long on the university pc...

scalinglawscutoff.mw

Ok, I can make the loops, but as seen on the worksheet, it doesn't return different values for different parameters. Also, for some reason this values are wrong. I did just for rc as you can see, but this happens to every quantity. For p0 and rho0 I don't even get numerical values. I interrupted the computation because it was taking too long on the university pc...

scalinglawscutoff.mw

Thanks I think I got now. I was just confused on how and where to include the loops. Thanks!

Thanks I think I got now. I was just confused on how and where to include the loops. Thanks!

@Markiyan Hirnyk Well, this was the old code, where I just needed to solve for one value for the parameters. What I'm trying to find out is an easy way to solve for several values of theta0, beta0 and thetaR without having to run dsolve every time. Then, somehow, I could call these solutions and extract stuff from then. 

I mean like a for loop

for beta from 1 to 100 do

 for theta0 from 1 to 100 do

    for thetaR from 1 to 100 do

       dsolve(system with initial conditions)

       plot(whatever)

       calculate all values I want 

    end do etc. etc.

then I'd have a file or a worksheet with the plots, and the calculations I want. But I'm not able to do that. The best thing I found is using the parameters option on dsolve and then using each value for the parameters I want, but this takes a lot of time since I have to calculate everything every time.

Hope this made things clearer