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]);