function [X,Y]=traject(func,z0,param,N,crit,maxit,epsilon);
%TRAJECT plots the trjectories of the following first-order system of 
%  differential equations:
%
%                    dx/dt=f(x,y,param)           (*)
%                    dy/dt=g(x,y,param)
%
%  INPUTS:
%    (i) func = a string giving the name of the m-file function with input
%               arguments z=[x;y] and param and output [f;g] (using notation
%               of (*) )
%   (ii) z0 = (2x1) starting point for trajectory (=[x0;y0])
%  (iii) param = parameters of 'func'.m
%   (iv) N = number of points to plot
%    (v) crit = criterion when recursively solving
%   (vi) maxit = maximum number of iterations for any recursive solving
%
%  OUTPUT: trajectories plotted
%
zbar=newton(func,z0,param,crit,maxit);
xbar=zbar(1); ybar=zbar(2);
a=xbar-epsilon(1);
for i=1:2*N+1;
  x(i)=a+(i-1)*epsilon(1)/N;
end;
for i=1:2*N+1;
  V0=ybar;
  W0=ybar;
  for j=1:maxit;
    eval(['[f,fj]=',func,'([x(i);V0],param);']);
    V1=V0-f(1)/fj(1,2);
    if max(abs(V1-V0))<crit; break; end;
    V0=V1;
  end;
  if j>=maxit; 
    sprintf('Iteration limit of %g reached when solving f(xbar,y(i))=0',maxit)
  end;
  for j=1:maxit;
    eval(['[f,fj]=',func,'([x(i);W0],param);']);
    W1=W0-f(2)/fj(2,2);
    if max(abs(W1-W0))<crit; break; end;
    W0=W1;
  end;
  if j>=maxit; 
    sprintf('Iteration limit of %g reached when solving g(xbar,y(i))=0',maxit)
  end;
  y1(i)=V1;
  y2(i)=W1;
end;
plot(x,y1,x,y2); pause;
for i=1:20;
  Z=input('Input starting point for trajectory (or [-99;-99] if done): ');
  X(1)=Z(1); Y(1)=Z(2);
  if X(1)==-99 & Y(1)==-99; break; end;
  for j=1:50;
    eval(['f=',func,'([X(j);Y(j)],param);']);
    X(j+1)=X(j)+f(1);
    Y(j+1)=Y(j)+f(2);
  end;
  plot(x,y1,x,y2,X,Y); pause;
end;

