solve({eq1,eq2},{d,s});

{d = 2*a*c^2+2*a^2*b+2*a^2*c+4*a*c*b+2*b^2*a+2*b^2*c+2*b*c^2-a-b-c, s = 2*a*c^2+2*a^2*b+2*a^2*c+4*a*c*b+2*b^2*a+2*b^2*c+2*b*c^2}

@philbertosaroni 

Let me try to clarify what the problem is, and please correct me where I misunderstand.

You have 6 non-linear differential equations in 6 variables. While the solution may not be unique, given fixed parameters and an appropriate number of boundary conditions, you should be able to simulate numerically a trajectory, even if only over a short period of time. It's possible, in general, that your system might go off to infinity or experience a discontinuity, etc.. but hopefully it won't for some of the parameter  values you are interersted in. 

You have set up 6 boundary conditions defined at the origin (t=0), i.e. initial conditions. You will need all 6 of them to simulate a solution. If you treat one of those as a parameter, then you can simulate a list of solutions for each value of the parameter you are willing to entertain. If you treat several initial conditions as parameters, you will have even more trajectories to look at. 

You do not need to "shoot" since you are specifying initial conditions and no final condition.

I suspect that what you mean to do, rather, is to simulate a saddle-path type system with a small number (unique?) of converging trajectories and an infinity of diverging ones, such that the initial conditions you have left unspecified must be selected to reach (if possible) the final conditions.  Is the ode system derived from an optimal control problem with state and control variables?

So, if my guess is correct, you'd need to specify where the system should end, perhaps a "stationary state" of the sytem (or a corner of the system, or whatever), so you would first solve numerically for those final conditions and use them to replace the unspecified initial conditions. You will need 6 conditions (x initial conditions and 6-x final conditions). You will also need to specify when to end the calculations, e.g. a long "infinite" time horizon or a "finite" range of time (t=0..10).

The boundary problem solvers work using shooting-like methods or other methods that you do not need to code yourself. You just feed the mixed boundary conditions and Maple will select a bc method.

What I have found works well with this kind of problem is this: Once you have your stationary state, where the system is supposed to end up asymptotically as time goes to infinity, then treat that as if it were your initial state, and then run the system backwards (negative time). With this approach you are solving an initial value problem, which seems to be easier to do: there are a wide range of numerical methods for "initial value problems", not so much flexibility for "boundary value problems". 

If you need a shooting algorithm, the Maple expert is Doug Meade, you can google your way to his shooting algorithms.

you could save the content as simple text, say as "input.txt" and then from within a worksheet:

 read "input.txt";

One approach is to get the solution of your system in the form of a list of points and then find the maximum from that list. The answer will be more precise the finer the grid you choose.

There may be other, more efficient methods. This is what I got:


#rewriting your input slightly. Note: I recommend you use 1D input in classic/standard worksheet.
restart; 
beta := 0.660707762617268e-4:
alpha := .229305575688023:
ode1 := diff(s(t), t) = -beta*s(t)*i(t):
ode2 := diff(i(t), t) = beta*s(t)*i(t)-alpha*i(t):
ode3 := diff(r(t), t) = alpha*i(t):
sys := [ode1,ode2,ode3];
ini := [i(0) = 50, r(0) = 0, s(0) = 10000];

sol := dsolve([op(sys),op(ini)], 'numeric', 'output' = listprocedure);

ps := plots:-odeplot(sol, [t, s(t)], t = 0 .. 50, 'colour' = red):
pi := plots:-odeplot(sol, [t, i(t)], t = 0 .. 50, 'colour' = green):
pr := plots:-odeplot(sol, [t, r(t)], t = 0 .. 50, 'colour' = blue):
plots:-display([ps,pi,pr]);

# eval(i(t),sol); plot(%, 0 ..50, 'discont'=true);# alternative way to plot

L := [ seq( eval([t,i(t)],sol) (k), k = 0 .. 50) ]:
max(L);

2902.84847845621

L := [ seq( eval([t,i(t)],sol) (k), k = 0 .. 50, 0.1) ]: # finer grid
max(L);

 2906.47960109601

1. main suggestion: use indentation to make sure every "if" has an "end if;" and every "do" has an "end do;"
2. don't forget to place an : or ; after each execution group
3. don't use protected symbols, like sum, product, add, etc..
4. no colon or semicolor after proc(), but instead after the locals are declared.
5. debug/test your proc by checking what each line does, this can be done by commenting out and using print or returning the current value

sum := 1;
Error, attempting to assign to `sum` which is protected

To illustrate indenting and debugging, consider:

myproc := proc()
local a,b,c,d,i,j,k,s,p,l,r,u;
r:=7.11: u:=0.1: # for debugging purposes, change the values here
for i by u to r do
     a := i;
     for j by u to r do
         b := j;
         for k by u to r do
             s := a+b+c+d;
             p := a*b*c*d;
             c := k;
             d := r-a-b-c;
             if p=s and s=r then
                 print(`When a=`,a,`b=`,b,`c=`,c,`and d=`,d);
                 print(`a+b+c+d=a*b*c*d=r`);
             else print("not working") # testing here
             end if;
         end do;
     end do;
end do;
end proc:

myproc();

Markiyan showed you how to get the results.

it depends what you mean by incrementing by 10, if you mean the exponent then that's exactly what Markiyan did with L. You can let the index run from -10 instead of 0. If you mean the numbers themselves, you'll run into memory issues, I think.

restart;
T := x -> 6.58*10^38*10^(-25)/x^3:
L := [seq(10^j, j = -10 .. 10)]:
# L := Array([seq(k, k = -10^7 .. 10^7), 10^2]); #error:  Error, object too large
Array(%id = 73572924)
plots:-pointplot( map(x -> [x, T(x)] , L)
, 'axis' = [ mode = log ]
, 'view' = min(L) .. max(L)
, 'symbol' = solidcircle
, 'symbolsize' = 12
, 'color' = red
) ;

# preliminary:
# seq(T, x in [10^10, 10^9, 10^8, 10^7, 10^6, 10^5, 10^4, 10^3, 10^2, 10, 10^0]);

#plotting the (x,T(x)) pairs  
plots:-pointplot( [seq([x,T], x in [10^10, 10^9, 10^8, 10^7, 10^6, 10^5, 10^4, 10^3, 10^2, 10, 10^0])]);

# controlling the view
plots:-pointplot( [seq([x,T], x in [10^10, 10^9, 10^8, 10^7, 10^6, 10^5, 10^4, 10^3, 10^2, 10, 10^0])]
     , 'view' = [default, 10^(-11) ..0]);

My contribution here is to give the perspective of an average user (student?) confronted with the problem. It appears to be very difficult to get the solution with Optimization:-Maximize, even tweaking the options, unless you pretty much know where the solution is, and even then ... 

Here are some steps that take you to the right place, but the journey is fraught with peril.

It's natural to guess that the solution might be on the 1..7 range, from there trial and error... it takes a few steps and some good intuition...

Note the erroneous solution:  [2.23606797749978981, [x = 1., y = 2.]]

restart;
f := sqrt((x-1)*(y-x))+sqrt((7-y)*(1-x))+sqrt((x-y)*(y-7)):

Optimization:-Maximize(f
, x = 1 .. 7
, y = 1 .. 7
, 'initialpoint' = [x=1,y=2]
) ;

Warning, undefined value encountered
Error, (in Optimization:-NLPSolve) number expected for float[8] parameter, got HFloat(HFloat(undefined))+HFloat(HFloat(undefined))*I


Optimization:-Maximize(f
, x = 2 .. 7
, y = 2 .. 7
, 'initialpoint' = [x=1,y=2]
) ;

Error, (in Optimization:-NLPSolve) complex value encountered


Optimization:-Maximize(f
, x = 1 .. 7
, y = 2 .. 7
, 'initialpoint' = [x=1,y=2]
) ;

Warning, no iterations performed as initial point satisfies first-order conditions
[2.23606797749978981, [x = 1., y = 2.]]


Optimization:-Maximize( subs(y=2,f)
, x = 1 .. 7
, y = 2 .. 7
, 'initialpoint' = [x=1,y=2]
) ;

Warning, no iterations performed as initial point satisfies first-order conditions
[2.23606797749978981, [x = 1., y = 2.]]


Optimization:-Maximize( subs(x=1,f)
, x = 1 .. 7
, y = 2 .. 7
, 'initialpoint' = [x=1,y=2]
) ;

[3., [x = 1., y = 4.]]

There's no way you could be confident to have found the max after all this.

But, as Axel pointed out, the plots do help.

If I'm not mistaken, the system analyzed by Preben (http://www.mapleprimes.com/view.aspx?sf=140842/448823/MaplePrimes12-11-2.mw) has a zero eigenvalue.

LinearAlgebra:-Eigenvectors(J(-1,y0,0));

Vector(3, {(1) = 0, (2) = -1, (3) = 6}), Matrix(3, 3, {(1, 1) = 0, (1, 2) = -1/y0, (1, 3) = -(1/7)/y0, (2, 1) = 1, (2, 2) = 1, (2, 3) = 17/21, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1})

[the bit that says (1, 1) = 0]

The system's stability can therefore be characterized by a center manifold, which could be center-stable or center-unstable, as the case may be. I'm not sure if the plot is enough to settle the matter. I'm not saying Preben's analysis isn't correct, I'm just saying that in some cases the plots can be deceptive and caution is needed.

And by the way a big thank you to Preben for posting this great worksheet.

Here is an analyis of dynamic systems that have a division by zero. It's rather subtle.

http://www.imati.cnr.it/~spinolo/BianchiniSpinolo4.pdf

I attach an adaptation of Preben's code to examples from the Bianchini & Spinolo (2011) article. The plots do not easily settle the matter. You can see that by playing around with the indeterminate critical ("initial") value.

You could apply their analysis to prove stability/instability for your system.

Incidentally, I was not able to export the animations as gif "programmatically" with plotsetup (but the diabolical mouse could), any idea what that is about? unsupported feature? bug? (I haven't tested on Maple 16.02)

I upload the worksheet, it's written on Maple 15.

In your second system escorpsy (I can no longer find it, did you remove it?) you have non-zero eigenvalues, so the usual Jacobian analysis is valid.

Bianchini-Spinolo-20.mw

Inserted animations: Are these guys stable?

you could declare your initial value as a parameter, this is explained here:

http://www.mapleprimes.com/posts/80873-Dsolve-With-Parameters

Not sure exactly what you mean. It depends on how you use solve. Say you have a quadratic with 2 solutions.

sol := [solve(x^2+x+1,x)];
sol[1];
sol[2];

As far as I know the order of the solutions is not necessarily consistent across worksheet executions, so you do get a list but the order of the list could change from one time to the next. 

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

I think you should cut-paste your unrelated question and create a new question for it, more efficient, more useful.

There are several things you can do, it may help to see a minimum working example of your worksheet to give advice. One approach is to save the data obtained from the simulations and to import that into your other worksheets.

If your worksheet is slow, it's quite possible that it could be optimized. Do you use the compile option? Show the basic layout of your dsolve commands if you want some advice on that.

One way to compute transformations of the original data is to incorporate them into an extended system and have dsolve deal with it right away, e.g. [x(t),y(t),z(t)], where z(t) is a transformed variable obtained from x(t) and y(t). This can be a better approach, depending on specifics.



Below are two procedures I use all the time.

exportData := proc(p::evaln)
local thedata, thedir, thename, i, n :
thedata := plottools:-getdata~([eval(p)]) :
thedir := cat(currentdir(),kernelopts(dirsep),data,kernelopts(dirsep)) :
n := nops([eval(p)]);
for i from 1 to n do
thename := cat(convert(p,string),"_",i,".txt") :
ExportMatrix(cat(thedir,thename), op([eval(i),-1],thedata), delimiter=",");
end do;
end proc:

importData := proc(p::evaln)
local thedir, thename, i, n :
thedir := cat(currentdir(),kernelopts(dirsep),data,kernelopts(dirsep)) :
n := nops([eval(p)]);
for i from 1 to n do
thename := cat(convert(p,string),"_",i,".txt") :
ImportMatrix(cat(thedir,thename), delimiter=",");
end do;
end proc:

If you have undefined data in your simulations (I often do), the following may help clean up and significantly reduce the size of the data/figures:

exportCleanData := proc(p::evaln)
local thedata, thedir, thename, M, U, i, j, m, n :
thedata := plottools:-getdata~([eval(p)]) :
thedir := cat(currentdir(),kernelopts(dirsep),data,kernelopts(dirsep)) :
U := NULL :
n := nops([eval(p)]) :
for i from 1 to n do
thename := cat(convert(p,string),"_",i,".txt"):
M := op([eval(i),-1],thedata):
m := LinearAlgebra:-RowDimension(M):
for j from 1 to m do
if has(LinearAlgebra:-Row(M,eval(j)), HFloat(undefined))
or has(LinearAlgebra:-Row(M,eval(j)), Float(undefined))
then U := U,j :
end if :
end do :
ExportMatrix(cat(thedir,thename), LinearAlgebra:-DeleteRow(M,[U]), delimiter=",") :
end do :
end proc :

one way is to use fprintf, here is how I've often exported:

fd := fopen(cat(theplace,thename,".m"), WRITE) :
fprintf( fd, %s, CodeGeneration[Matlab](thedata,'output'=string) ):
fclose(fd) :

Check CodeGeneration for other options and fprintf for other formats.

 (If I understand your question correctly), you could create a Matrix or Array of the data you want to export and then do something along these lines:



If N denotes your N'th simulation, say, N:=10, with a folder for each simulation (my practice):

N:=10:
theplace := cat(currentdir(),kernelopts(dirsep),convert(N,string),kernelopts(dirsep)); # name of the directory path
FileTools:-MakeDirectory(theplace); # create directory if it doesn't exist

thedata := Matrix([1,2,3]); #create data Matrix here

ExportMatrix(cat(theplace,"thedata",convert(N,string),".m")
, thedata, delimiter="&", format=rectangular, mode=ascii):


Check ExportMatrix for options, you can save in Matlab format. Here I added a delimiter &, which is convenient if you're going to copy-paste into a LaTeX table, just an example.