:

## second follow-up on my second blog post

Maple

Below is the latest version of the code to draw iso-chrone lines and a salvo of arrows onto the phase diagram of a two-dimensional system of ordinary differential equations.

I have greatly benefited from inputs by Robert Israel (who wrote the first incarnation of the procedure), Joe Riel, and pagan. A big thankyou!

The procedure is sufficiently developed for my current purpose, so I don't plan to modify it much in the near future.

Tested on Maple 13/ Classic. The last plot combines the isochrones and the salvo.

```
restart;

# Using the undocumented _Env_in_maplet variable
getarrows0 := proc(x0, y0, T, x1, y1)
local DEsol,tm,j,wl;
global x,y,t,integ;
_Env_in_maplet := true;
integ('initial' = [x0,y0]);  # reinitialize
integ('parameters'=[x1,y1]); # assign parameters
tm:= floor(subs(integ(T),t));
seq(subs(integ(j), [t, x(t), y(t)]), j = 0 .. tm);
end proc:
```
```# Not using  _Env_in_maplet and less efficient
getarrows := proc(x0, y0, T, x1, y1)
local DEsol,tm,j,wl;
global x,y,t,integ;
wl := interface('warnlevel'=0); # silence event warnings during integration
try
integ('initial' = [x0,y0]);  # reinitialize
integ('parameters'=[x1,y1]); # assign parameters
tm:= floor(subs(integ(T),t));
finally
interface('warnlevel'=wl); # restore warnlevel even if an error occurs
end try;
seq(subs(integ(j), [x(t), y(t)]), j = 0 .. tm);
end proc:

xdot := diff(x(t),t) = x(t)-y(t):
ydot := diff(y(t),t) = y(t)*(1-x(t)):
integ := dsolve({xdot,ydot,x(0)=0,y(0)=0}
, numeric
, 'events' = [[y(t)-y1,'halt'],[x(t)-x1,'halt'],[x(t),'halt']]
, 'output=procedurelist'
, 'parameters'=[x1,y1]
):
```
```T := 10:
xmax := 3:
ymax := 3:
DEsol:= dsolve({xdot,ydot,x(0)=1-1e-10,y(0)=1-1e-10},numeric,output=listprocedure):
t1 := fsolve(rhs(DEsol(t)[2])=0.05, t=-100..0):
x0 := rhs(DEsol(t1)[2]):
y0 := rhs(DEsol(t1)[3]):

st := time[real]():
y0a := [ seq( seq( y0*(1+i/10*(1e-1)^j),i=-90..99 ), j=1..10) ]:
L := map2(getarrows0, x0, y0a, T, xmax, ymax):
L2 := seq( select( s->s[1]=t, L ), t=0..T ):
iso := t-{op(L2[t+1,..][..,2..3])}:
init1a := L[..,2..3]:
time[real]() - st;

p1 := seq( plots:-pointplot( iso(t), connect=true, thickness=3 ),t=0..T ):
```
```plots:-display(p1);

p2 :=
DEtools[DEplot](
{xdot,ydot},
{x(t),y(t)},
t=-20..20,
[[x(0)=0.9999,y(0)=0.9999],[x(0)=1.0001,y(0)=1.0001]],
scene=[x(t),y(t)],
x=0..xmax, y=0..ymax,
arrows=smalltwo,
dirfield=init1a,
linecolor=black,
thickness=2,
numpoints=1000
):

plots:-display(p2, view=[0..2,0..2]);

p3 :=
DEtools[DEplot](
{xdot,ydot},
{x(t),y(t)},
t=-20..20,
[[x(0)=0.9999,y(0)=0.9999],[x(0)=1.0001,y(0)=1.0001]],
scene=[x(t),y(t)],
x=0..xmax, y=0..ymax,
arrows=none,
linecolor=blue,
thickness=5,
numpoints=1000
):

plots:-display(p3):
```
```y0b := [ seq( seq( y0*(1+i*(1e-1)^j),i=-3..3 ), j=1..5) ]:
Lb := map2(getarrows0, x0, y0b, T, xmax, ymax):
init1b := Lb[..,2..3]:

p4 :=
DEtools[DEplot](
{xdot,ydot},
{x(t),y(t)},
t=-20..20,
[[x(0)=0.9999,y(0)=0.9999],[x(0)=1.0001,y(0)=1.0001]],
scene=[x(t),y(t)],
x=0..xmax, y=0..ymax,
arrows=smalltwo,
dirfield=init1b,
color=red,
linecolor=blue,
thickness=5,
numpoints=1000
):

plots:-display(p4):
```
```plots:-display(p1,p4, view=[0..2,0..2]);
```

﻿