@PatrickT 

This is the way acer's code deals with initial conditions: (N.B. I'm copy-pasting-editing so some risk of typos, but hopefully the gist will get across)

  h := 7; # maximum initial values of x1
  s := h/10; # size of step
  k := (-2)*(1/(h-1)); 
  d := abs(k)+1;
  INI := ii -> [x[1](0) = ii, x[2](0) = k*ii+d]: INI(ii);

It then uses the option parameters in dsolve to handle the parameter ii, e.g.

ZZ := dsolve([SYS, INI(ii)], [x[1](t), x[2](t)], 'parameters' = [ii], 'numeric');

and then, after defining the "master" procedure, it uses seq over the parameter ii (as well as other parameters I ignore here), zooming on the relevant piece of code:

master := subs(dummy = eval(ZZ), proc (II); dummy('parameters' = ['ii' = II]); dummy; end proc);
seq(plots:-odeplot(master(i), [x[1](t), x[2](t)], t=-1 .. 1), i = 1 .. h, s)

What if instead of sequencing over ii, you want to feed in a predefined list of ii. The adaptation is quite straightforward. Say you have (as I did):

Then proceed as before with

i = candidates

instead of

i = 1 .. h, s

But then, what if you want to control the color scheme? The following worked. I first create a list of colors I like:

colorlist := [black, red, blue, yellow]: # if you have a fixed number of "candidates" or if you have a flexible number of "candidates" N:
colorlist := [seq(COLOR(HUE,k/N),k=1..N)]:                        # Color Sequence for Plots
colormap := i -> colorlist[ListTools:-Search(i,candidates)]:            # Map Candidate to Color

and then add the option,

'colour' = colormap(i)

where the index i is the same as the one in master(i)

so to sum up, you have something like:

master := subs(dummy = eval(ZZ), proc (II); dummy('parameters' = ['ii' = II]); dummy; end proc);
p := seq(plots:-odeplot(master(i), [x[1](t), x[2](t)], t=-1 .. 1)
  , 'colour' = colormap(i)
)
  , i = candidates
):
plots:-display(p);

@acer 

This is a great tip acer, thanks. I notice tremendous efficiency improvements.

I have used the option remember in the past, with great effect, but somehow did not make the connection here with the output of dsolve. Even if I had though, I would have had a hard time implementing it.

subsop(3=remember,eval(x(t), SOL)):

Not exactly your user friendly piece of code! Google returns only a very small instances of its use. After fooling around a little I understand what it does. First, eval(x(t), SOL); returns a procedure. And as such it can take options. As it happens, the options are the third argument, which can be seen by examining the output of op(eval(x(t), SOL));

Is it always the case that the third argument of a procedure is the "options"? does

subsop(3=something

writes over existing options? or does it merely add extra options? (this is a question to myself, I can easily experiment to find out).

Is there a more intuitive way to access and modify a procedure's options? like ModifyProc(myproc, 'option'=remember);

So thanks again acer for this tip. I will now try

subsop(3=remember, mybrain);

This is just to say, in response to an email from the site admin, that yes it is okay to branch this sub-discussion off. And this comment could be removed in the process. With thanks. Patrick.

So, if I understand correctly, posts that get voted up are listed higher up on the thread. This is a very bad idea, in my opinion. Here's the sort of thing that happens: an answer appears above a question. As a result, it gets very difficult to follow the conversation.

The answer got a vote up, while the question didn't, as you would expect. Result: the answer appears before the question.

@acer 

The other thing I need to do with the output is plot functions of the solution, as opposed to the solution itself. If x(t) and y(t) are the variables, I need to plot things like (x(t))^2 or 1/y(t) or more complicated functions of these.

I used to do it this way:

SOL := dsolve(blabla):

xt := eval( x(t), SOL ):
yt := eval( y(t), SOL ):
zt := eval( z(t), SOL ):
ft := eval( f( xt, yt, zt ) ):

plots[odeplot]( SOL, [ ft(t), xt(t) ] ):   # x versus f(x,y,z)

where f is a known function.

But this is not terribly efficient if I want to plot several plots like the one above, because every call to plots[odeplot]( SOL,...) reruns dsolve.

@acer 

I have been able to adapt your code to the problem that I'm currently working on, a 3-D system for which I simulate the dynamics on a trajectory converging on a stationary point for a range of parameter values (the stationary point is characterized by a center manifold and somewhat stiff). Having the code written in your way is convenient, it saves me from having to rewrite and rerun stuff. I cannot really comment on performance, I trust your judgement on that.

I have been able to extend your code to dfieldplot3d, without any difficulty. This would be particularly useful if I could project the 3D vector field onto a 2D plane or surface that I would specify. I have not experimented on this yet.

My main use so far has been to plot a list of trajectories of the solution of a 3D system but projected onto 2D and without the vector field, e.g. against time  [t,x(t)] or against each other [x(t),y(t)].

There is a limitation I have run into from using 'parameters' in dsolve. Perhaps there is a workaround.

I sometimes like to output of dsolve to be projected onto an array of points, usually like this:

T:=100: N:=20: pts := Array(1..N, i->(i-1)*abs(T)/N):
dsolve( [SYS, INI], VAR, numeric, output = pts);

Unfortunately, neither 'events' nor 'parameters' can be used with array-style output.

I can use plots:-odeplot with the option "style=point" but then I do not know how to control the spacing and number of points of the output.

I'll report back if I have further questions / answers about the above. Thanks again acer, you've been very helpful!

To accommodate 3-D systems as well as 2-D systems, the code may be rewritten as follows. At this time, one needs to select the number of variables N (either 2 or 3) and modify the vector field f accordingly, and also toggle dfieldplot to dfieldplot3d (that could be made automatic within a procedure, but at this time the bother does not seem worthwhile):

N:=3: # number of variables

# 3-D Vector field, to illustrate:
f := proc (x,a)
  options operator, arrow;
  [
     -(1-x[1]^2)*(2*x[2]^2+(1/2+(3/2)*a[2])*(1-x[2]^2)-x[1]*x[2])
   , -(2-2*x[2]^2-(1/2+(3/2)*a[2])*(1-x[2]^2))*x[1]*x[2]-(1-x[1]^2)*(1-x[2]^2)-(3/2-(3/2)*x[2]^2)*(1-a[2])*a[1]
   , x[1]+x[2]
  ]:
end proc:
F := f(x,a);

... insert same code as posted earlier...

# Options
a2:=-.2:
xrange := [1 .. 5, -1.1 .. 1.1, -1 .. 1]:
theview := [1 .. 6, -1 .. 1, 0 .. 1]:
A1range := -4 .. 10:
thegrid := [10, 10, 10]:

# Create animation frames
allframes :=
[
  seq(
    plots:-display(
      seq(plots:-odeplot(master(A1, a2, i), VAR, t=-1 .. 1, 'numpoints' = 500, 'thickness' = 3), i = 1 .. h, s),
      plots:-fieldplot3d(                     #2-D:fieldplot, 3-D:fieldplot3d
        eval(f(x,A),[A[1]=A1,A[2]=a2])
         , seq(x[i]=xrange[i],i=1..N)
         , 'fieldstrength' = 'fixed'          #maximal,average,fixed,log,scaleto(v)
         , 'grid' = [seq(thegrid[i],i=1..N)]  #default=[20,20]
        , 'arrows' = SLIM                    #LINE,THIN,SLIM,THICK,`3-D` # CAPITAL LETTERS for fieldplot3d
      )
    , view = [seq(theview[i],i=1..N)]
    , 'axes' = 'box'
  ), A1 = A1range )
]:

I'm sure I need your code acer. I have some pretty inefficient simulations that would greatly benefit from using your approach. So I've tried to understand the general gist of the code. My objective was to write a template that I could easily adapt to my own systems, without having to fiddle too much at each adaptation.

I have this system of 3 variables, can your code be adapted to a 3D system?

With a view to adapting the code to my own needs, I have changed the notation a little. This was mostly an exercise to get to learn your code. The next step perhaps would be to write the code with Vectors or Lists, so that the number of parameters, in particular, could be flexible.

As a memento to myself if nothing else, here's the adaptation, which produces the same thing as the first incarnation of your code (apart from minor changes):

N:=2: # number of variables
M:=2: # number of parameters

# Vector Field
f := proc (x,a)
  options operator, arrow;
  # Option arrow, in conjunction with option operator, indicates that the operator was initially entered using the -> notation.
  # The use of option arrow also disables Maple simplification rules that add any non-local names in the operator expression to the parameter list.
  # This option has no meaning for modules.
  [
   -(1-x[1]^2)*(2*x[2]^2+(1/2+(3/2)*a[2])*(1-x[2]^2)-x[1]*x[2]),
   -(2-2*x[2]^2-(1/2+(3/2)*a[2])*(1-x[2]^2))*x[1]*x[2]-(1-x[1]^2)*(1-x[2]^2)-(3/2-(3/2)*x[2]^2)*(1-a[2])*a[1]
  ]:
end proc:
F := f(x,a);

# Dynamic System
sys := (x,a,t) -> [seq( diff(x[j](t),t) = eval(f(x,a),[seq(x[i]=x[i](t),i=1..N)])[j] , j=1..N )]:
SYS := op(sys(x,a,t));

# Variables
var := (x,t) -> [seq(x[i](t),i=1..N)]:
VAR := var(x,t);

# Set initial values for x1 and x2
  h := 7; # maximum initial values of x1
  s := h/10; # size of step
  k := (-2)*(1/(h-1));
  d := abs(k)+1;

# Initial Values
ini := (x,j) -> [x[1](0) = j, seq(x[i](0) = k*j+d,i=2..N)]:
INI := op(ini(x,ii)): INI;

# Parameters
par := a -> [seq(a[i],i=1..M),ii]:
PAR := par(a);

# Master
ZZ := dsolve([SYS,INI], VAR, 'parameters' = PAR, 'numeric');
master := subs(dummy = eval(ZZ),
  proc (A1, A2, II)
    dummy('parameters' = ['a[1]' = A1, 'a[2]' = A2, 'ii' = II]);
    dummy
  end proc);

# Options
a2:=-.2:
x1range := 1 .. 15:
x2range := -1.1 .. 1.1:
theview := [1 .. 13, -1 .. 1]:
A1range := -4 .. 10:
thegrid := [20, 20]:

# Create animation frames
allframes :=
[
  seq(
    plots:-display(
      seq(plots:-odeplot(master(A1, a2, i), VAR, t=-1 .. 1, 'numpoints' = 1000, 'thickness' = 3), i = 1 .. h, s),
      plots:-fieldplot(
        eval(f(x,A),[A[1]=A1,A[2]=a2])
        , x[1] = x1range
        , x[2] = x2range
        , 'fieldstrength' = 'fixed'  #maximal,average,fixed,log,scaleto(v)
        , 'grid' = thegrid           #default=[20,20]
        , 'arrows' = slim            #line,thin,slim,thick
      )
    , view = theview
    , 'axes' = 'box'
  ), A1 = A1range )
]:

plots:-display(allframes, 'insequence' = true, color = blue);

Another respect in which the site's potential is not optimized is that often stuff uploaded on the site in the past, e.g. pictures, is no longer available, an example I just came across, from December 2008:

http://www.mapleprimes.com/questions/38142-How-To-Plot-The-Function-Dsolve-Function#comment68262
http://www.mapleprimes.com/files/97_Picture%204.png

When I clicked on it today, the picture (supposedly uploaded to the site) was no longer available.

nice one Alec! (good to see you're actively posting again)

nice one Alec! (good to see you're actively posting again)

Great stuff acer, I'm often looking for ways to improve the efficiency of my ODE simulations, will put this tutorial to great use.

As a side remark, if the (not needed) plottools package happens to be loaded before running the code, it will interfere and cause an error, as I have experienced.

thanks pagan,

so the syntax is events=list of list with the trigger and the event separated by a comma.

This is actually quite intuitive.

Side remark:

I misinterpreted the syntax because the Help example contains two triggers and hence a nested list, a simplified version of the example from the help file:

events = [[[y(t), 0 < diff(y(t), t)], halt]]

thanks pagan,

so the syntax is events=list of list with the trigger and the event separated by a comma.

This is actually quite intuitive.

Side remark:

I misinterpreted the syntax because the Help example contains two triggers and hence a nested list, a simplified version of the example from the help file:

events = [[[y(t), 0 < diff(y(t), t)], halt]]

There seems to be general agreement that certain things could be improved. I wouldn't go so far as to say that the site has become "unusable," that's perhaps a bit of an exageration (understandable under the circumstances). I think (I could be wrong, my memory's awful) spam started under the old mapleprimes (is that right?), I seem to remember something to the effect that the new mapleprimes was going to fix the spam issue, well it didn't fix it, but I don't think it's fair to blame spam on it either. Spam's a big problem all over the net.

I think we can all agree that it would be great if 1) search could be improved (it should at least match google's ability to catch stuff within mapleprimes), if 2) messages could be ordered on the fly according to date (especially), users, etc.with some filters, and for this the "comments" should not be treated as some sub-standard message that doesn't get listed if 3) the voting and ranking system could be simplified and clarified, for instance by eliminating the down-voting that upsets nearly everyone, by allowing comments to be voted on, by clarifying the way the counters are incremented, if 4) the email notification could be improved (Alejandro, I use it too, but at times I get 5 or 6 notifications for the same update and sometimes none, not as "useable" as it should be, I agree). And a few other things too that I can't remember right now.