function output = pplane5(action,input1,input2,input3)

% PPLANE5  is an interactive tool for studying planar autonomous systems of
%	differential equations.  When PPLANE5 is executed, a PPLANE5 Setup
%	window is opened.  The user may enter the differential
%	equation and specify a display window using the interactive
%	controls in the Setup window.  Up to 4 parameters may also be 
%	specified.  In addition the user is give a choice of the type of 
%	field displayed and the number of field points.
%
%	When the Proceed button is pressed on the Setup window, the PPLANE5
%	Display window is opened, and a field of the type requested is 
%	displayed for the system.  When the mouse button is depressed in 
%	the PPLANE5 Display window, the solution to the system with that
%	initial condition is calculated and plotted.
%
%	Other options are available in the menus.  These are
%	fairly self explanatory.  
%
%	This is version 5.0b10, and will only run on version 5 of MATLAB.

%  Copywright (c)  1995, 1997, 1998 John C. Polking, Rice University
%  Last Modified: March 13, 1997

startstr = 'pplane5';
if nargin < 1
  action ='initialize';
end

if get(0,'screendepth') == 1
   style = 'bw';
end

if strcmp(action,'initialize')
   
   % First we make sure that there is no other copy of PPLANE5
   % running, since this causes problems.	
   pph = findobj('tag','PPLANE5');
   if ~isempty(pph);
      qstring = {'There are some PPLANE5 figures open although they may be invisible.  ';...
            'What do you want to do?'};
      tstring = 'Only one copy of PPLANE5 can be open at one time.';
      b1str = 'Restart PPLANE5.';
      b2str = 'Just close those figures.';
      b3str = 'Do nothing.';
      answer = questdlg(qstring,tstring,b1str,b2str,b3str,b1str);
      if strcmp(answer,b1str)
         delete(pph);
         %pplane5;return
         eval(startstr);return
      elseif strcmp(answer,b2str)
         delete(pph);return
      else
         return
      end 
   end  % if ~isempty(pph);
   
   % Next we look for old files created by PPLANE5.
   
   oldfiles = dir('pptp*.m');
   kk = zeros(0,1);
   for k = 1:length(oldfiles)
      fn = oldfiles(k).name;
      fid = fopen(fn,'r');
      ll = fgetl(fid);
      ll = fgetl(fid);
      ll = fgetl(fid);
      fclose(fid);
      if strcmp(ll,'%% Created by PPLANE5')
         delete(fn)
      else
         kk = [kk;k];
      end
   end
   if ~isempty(kk)
      fprintf('The following files\n');
      for j = 1:length(kk)
         fprintf(['     ',oldfiles(kk(j)).name,'\n']);
      end
      fprintf('may have been created by an older version of PPLANE.\n');
      fprintf('If so they should be deleted.\n\n');
   end
   
   style = 'white';
   ppdir = pwd;
   ssize = get(0,'defaultaxesfontsize');
   npts = 20;
   solver = 'Dormand Prince';
   tol = 5e-4;
   stepsize = 0.1;
   if exist('ppstart','file')
     H = ppstart;
     if ~isempty(H)
       if isfield(H,'style')
	 style = H.style;
       end
       if isfield(H,'size')
	 ssize = H.size;
       end
       if isfield(H,'npts')
	 npts = H.npts;
       end
      if isfield(H,'solver')
	solver = H.solver;  
      end
      if isfield(H,'ppdir')
	ppdir = H.ppdir;  
      end
      if isfield(H,'stepsize')
	stepsize = H.stepsize;  
      end
      if isfield(H,'tolerance')
	tol = H.tolerance;  
      end
    end
  end
  if get(0,'screendepth') == 1
    style = 'bw';
  end
  
  ud.ssize = ssize;
  ud.ppdir = ppdir;
  comp = computer;
  if strcmp(comp,'PCWIN')
    ud.fontsize = 0.8*ud.ssize;
  else
    ud.fontsize = ud.ssize;
  end
   
  % Set up for the menu of systems.
  
  system.name = 'default system';
  system.xvar = 'x';
  system.yvar = 'y';
  system.xder = ' 2*x - y + 3*(x^2-y^2) + 2*x*y';
  system.yder = ' x - 3*y - 3*(x^2-y^2) + 3*x*y';
  system.pname1 = '';
  system.pname2 = '';
  system.pname3 = '';
  system.pname4 = '';
  system.pval1 = '';
  system.pval2 = '';
  system.pval3 = '';
  system.pval4 = '';
  system.fieldtype = 'arrows';
   system.npts = npts;
   system.wind = [-2 4 -4 2];
   
   system(2).name = 'linear system';
   system(2).xvar = 'x';
   system(2).yvar = 'y';
   system(2).xder = ' A*x + B*y';
   system(2).yder = ' C*x + D*y';
   system(2).pname1 = 'A';
   system(2).pname2 = 'B';
   system(2).pname3 = 'C';
   system(2).pname4 = 'D';
   system(2).pval1 = '2';
   system(2).pval2 = '2';
   system(2).pval3 = '-2';
   system(2).pval4 = '-3';
   system(2).fieldtype = 'arrows';
   system(2).npts = npts;
   system(2).wind = [-5 5 -5 5];
   
   system(3).name = 'vibrating spring';
   system(3).xvar = 'x';
   system(3).yvar = 'v';
   system(3).xder = ' v';
   system(3).yder = ' -(k*x + d*v)/m';
   system(3).pname1 = 'k';
   system(3).pname2 = 'm';
   system(3).pname3 = 'd';
   system(3).pname4 = '';
   system(3).pval1 = '3';
   system(3).pval2 = '1';
   system(3).pval3 = '0';
   system(3).pval4 = '';
   system(3).fieldtype = 'arrows';
   system(3).npts = npts;
   system(3).wind = [-5 5 -5 5];
   
   system(4).name = 'pendulum';
   system(4).xvar = '\theta';
   system(4).yvar = '\omega';
   system(4).xder = ' \omega';
   system(4).yder = ' -sin(\theta) - D*\omega';
   system(4).pname1 = 'D';
   system(4).pname2 = '';
   system(4).pname3 = '';
   system(4).pname4 = '';
   system(4).pval1 = '0';
   system(4).pval2 = '';
   system(4).pval3 = '';
   system(4).pval4 = '';
   system(4).fieldtype = 'arrows';
   system(4).npts = npts;
   system(4).wind = [-10 10 -4 4];
   
   system(5).name = 'predator prey';
   system(5).xvar = 'prey';
   system(5).yvar = 'predator';
   system(5).xder = ' (A - B*predator)*prey';
   system(5).yder = ' (D*prey - C)*predator';
   system(5).pname1 = 'A';
   system(5).pname2 = 'B';
   system(5).pname3 = 'C';
   system(5).pname4 = 'D';
   system(5).pval1 = '0.4';
   system(5).pval2 = '0.01';
   system(5).pval3 = '0.3';
   system(5).pval4 = '0.005';
   system(5).fieldtype = 'arrows';
   system(5).npts = npts;
   system(5).wind = [0 120 0 80];
   
   system(6).name = 'competing species';
   system(6).xvar = 'x';
   system(6).yvar = 'y';
   system(6).xder = ' (1 - x - y)*x';
   system(6).yder = ' (A - B*x - C*y)*y';
   system(6).pname1 = 'A';
   system(6).pname2 = 'B';
   system(6).pname3 = 'C';
   system(6).pname4 = '';
   system(6).pval1 = '4';
   system(6).pval2 = '2';
   system(6).pval3 = '7';
   system(6).pval4 = '';
   system(6).fieldtype = 'arrows';
   system(6).npts = npts;
   system(6).wind = [0 1 0 1];
   
   system(7).name = 'van der Pol''s equation';
   system(7).xvar = 'x';
   system(7).yvar = 'y';
   system(7).xder = ' M*x - y - x^3';
   system(7).yder = ' x';
   system(7).pname1 = 'M';
   system(7).pname2 = '';
   system(7).pname3 = '';
   system(7).pname4 = '';
   system(7).pval1 = '2';
   system(7).pval2 = '';
   system(7).pval3 = '';
   system(7).pval4 = '';
   system(7).fieldtype = 'arrows';
   system(7).npts = npts;
   system(7).wind = [-5 5 -5 5];
   
   system(8).name = 'Duffing''s equation';
   system(8).xvar = 'x';
   system(8).yvar = 'y';
   system(8).xder = ' y';
   system(8).yder = ' -(k*x + c*y + l*x^3)/m';
   system(8).pname1 = 'k';
   system(8).pname2 = 'c';
   system(8).pname3 = 'm';
   system(8).pname4 = 'l';
   system(8).pval1 = '-1';
   system(8).pval2 = '0';
   system(8).pval3 = '1';
   system(8).pval4 = '1';
   system(8).fieldtype = 'arrows';
   system(8).npts = npts;
   system(8).wind = [-3 3 -3 3];
   
   ud.c = system(1);	% Changed values.
   ud.o = system(1);	% Original values.
   ud.h = system(1);	% This will be the handles in the
   % setup window.
   
   ud.style = style;
  ud.npts = npts;  
   ud.settings.magn = 1.25;
   ud.settings.refine = 8;
   ud.settings.tol = tol;
   ud.settings.solver = solver;
  ud.settings.stepsize = stepsize;
   ud.system = system;
   
   switch style
   case 'black'
      color.temp = [1 0 0]; % red for temporary orbits
      color.orb = [1 1 0];  % yellow for orbits
      color.eqpt = [1 0 0];  % red for eq. pts.
      %    color.arrows = [0 1 1]; % cyan for arrows
      % color.arrows = [.5 .5 .9];  % purple for arrows
      color.arrows = .5*[1 1 1];  % gray for arrows
      color.narrows = .7*[1 1 1];  % gray for nullcline arrows
      color.tx = [1 1 0]; % yellow for xt plots & 3D plots
      color.ty = [1 0 0]; % red for yt plots
      color.ta = [1 0 0]; % red for axis plots
      color.sep = [.2,1,0]; % green for separatrices
      color.xcline = [1 1 1]; % white for xclines
      color.ycline = [1 0 1]; % magenta for yclines
      
   case 'white'
      color.temp = [1 0 0]; % red for temporary orbits
      color.orb = [0 0 1];  % blue for orbits
      color.eqpt = [1 0 0];  % red for eq. pts.
      color.arrows = 0.7*[1 1 1]; % gray for arrows
      color.narrows = .4*[1 1 1];  % gray for nullcline arrows
      color.tx = [0 0 1]; % blue for xt plots & 3D plots
      color.ty = [1 0 0]; % red for yt plots
      color.ta = [1 0 0]; % red for axis plots
      color.sep = [0,1,0];% green for separatrices
      color.xcline = [1 0 .75]; % magenta for xclines
      color.ycline = [1 .5 0]; %  orange for yclines
      
    case 'test'
      color.temp = [1 0 0]; % red for temporary orbits
      color.orb = [0 0 1];  % blue for orbits
      color.arrows = .7*[1 1 1]; % gray for arrows
      color.eqpt = [1 0 0];  % red for eq. pts.
      color.narrows = .4*[1 1 1];  % gray for nullcline arrows
      color.tx = [0 0 1]; % blue for xt plots & 3D plots
      color.ty = [1 0 0]; % red for yt plots
      color.ta = [1 0 0]; % red for axis plots
      color.sep = [0,1,0];% green for separatrices
      color.xcline = [1 0 .75]; % magenta for xclines
      color.ycline = [1 .5 0]; %  orange for yclines

    case 'display'
      color.temp = [1 0 0]; % red for temporary orbits
      color.orb = [0 0 1];  % blue for orbits
      color.arrows = .4*[1 1 1]; % gray for arrows
      color.eqpt = [1 0 0];  % red for eq. pts.
      color.narrows = .4*[1 1 1];  % gray for nullcline arrows
      color.tx = [0 0 1]; % blue for xt plots & 3D plots
      color.ty = [1 0 0]; % red for yt plots
      color.ta = [1 0 0]; % red for axis plots
      color.sep = [0,1,0];% green for separatrices
      color.xcline = [1 0 .75]; % magenta for xclines
      color.ycline = [1 .5 0]; %  orange for yclines
      
   case 'bw'
      color.temp = [1 1 1]; % white for everything
      color.orb = [1 1 1];  
      color.eqpt = [1 1 1]; 
      color.arrows = [1 1 1];
      color.narrows = [1 1 1];  
      color.tx = [1 1 1]; 
      color.ty = [1 1 1]; 
      color.ta = [1 1 1]; 
      color.sep = [1 1 1];
      color.xcline = [1 1 1];
      color.ycline = [1 1 1];
   end
   ud.color = color;
   
   ppset = figure('name','PPLANE5 Setup','numb','off',...
      'tag','PPLANE5','visible','off',...
      'user',ud);
   
   pplane5('figdefault',ppset);
   frame(1) = uicontrol('style','frame','visible','off');
   eq(1)=uicontrol('style','text',...
      'horizon','center',...
      'string','The differential equations.','visible','off');
   ext = get(eq(1),'extent');
   rr=ext(4)/10;
   
   texth =ext(4)+4;    % 19;   		% Height of text boxes.
   varw = 40*rr;		% Length of variable boxes.
   equalw =13*rr;		% Length of equals.(30)
   eqlength = 230*rr;		% Length of right hand sides of equations.
   winstrlen = 120*rr;	% Length of string boxes in display frame.
   left = 0;		% Left margin of the frames.
   frsep = 0;    %3;		% Separation between frames.
   separation = texth;	% Separation between boxes.
   
   dfigwidth =2*left + varw+equalw+eqlength+10;	% Width of the figure.
   dfigurebot = 30;	% Bottom of the figure.
   buttw = dfigwidth/3;
   qwind = [0,frsep,buttw,separation];	% Quit button
   rwind = [buttw,frsep,buttw,separation];	% Revert "
   pwind = [2*buttw,frsep,buttw,separation];	% Proceed "
   
   disfrbot = 2*frsep + separation;	% Display frame.
   disfrw = winstrlen + varw +10;
   disfrht = 5*separation + 10;
   disfrwind = [left, disfrbot, disfrw, disfrht];
   
   pfrbot = disfrbot + disfrht +frsep;	% Parameter frame.
   pfrw = dfigwidth -2*left;
   pfrht = 2*separation + 10;
   pfrwind = [left, pfrbot, pfrw, pfrht];
   
   defrbot = pfrbot + pfrht + frsep;	% Equation frame.
   defrw = pfrw;
   defrht = 3*separation + 10;
   defrwind = [left, defrbot, defrw, defrht];
   
   ffrbot = disfrbot;			% Field frame.
   ffrleft = left + disfrw + frsep;
   ffrw = dfigwidth -left - ffrleft;
   ffrht = disfrht; 
   ffrwind = [ffrleft, ffrbot, ffrw, ffrht];
   
   dfigureheight = defrbot + defrht +frsep;	% Height of the figure.
   
   set(ppset,'pos',[30 dfigurebot dfigwidth dfigureheight]);
   
   set(frame(1),'pos',defrwind);
   
   
   xname=[
      'ud = get(gcf,''user'');'...
         'Xname=get(ud.h.xvar,''string'');'...
         'set(ud.h.twind(1),''string'',[''The minimum value of '',Xname,'' = '']);'...
         'set(ud.h.twind(2),''string'',[''The maximum value of '',Xname,'' = '']);'...
         'ud.c.xvar = Xname;'...
         'ud.flag = 0;'...
         'ud.c.name = '''';',...
         'set(gcf,''user'',ud);'];
   
   yname=[
      'ud = get(gcf,''user'');'...
         'Yname=get(ud.h.yvar,''string'');'...
         'set(ud.h.twind(3),''string'',[''The minimum value of '',Yname,'' = '']);'...
         'set(ud.h.twind(4),''string'',[''The maximum value of '',Yname,'' = '']);'...
         'ud.c.yvar = Yname;'...
         'ud.flag = 0;'...
         'ud.c.name = '''';',...
         'set(gcf,''user'',ud);'];
   
   xder =[
      'ud = get(gcf,''user'');'...
         'ud.c.xder = get(ud.h.xder,''string'');'...
         'ud.flag = 0;'...
         'ud.c.name = '''';',...
         'set(gcf,''user'',ud);'];
   
   yder =[
      'ud = get(gcf,''user'');'...
         'ud.c.yder = get(ud.h.yder,''string'');'...
         'ud.flag = 0;'...
         'ud.c.name = '''';',...
         'set(gcf,''user'',ud);'];
   
   equationbot = defrbot + 5;
   eqlabelbot = equationbot + 2*separation;
   xbot = equationbot + separation;		% Bottom of x equation.
   ybot = equationbot;				% Bottom of y equation.
   
   lablen =200*rr;
   
   eqlableft = (dfigwidth-lablen)/2;
   eqleft = left + 5;
   
   set(eq(1),'pos',[eqlableft eqlabelbot lablen texth]);
   
   ud.h.xvar=uicontrol('pos',[eqleft, xbot, varw, texth],...
      'style','edit',...
      'horizon','right',...
      'string',ud.o.xvar,...
      'call',xname,...
      'visible','off');
   
   eq(2) = uicontrol('style','text',...
      'pos',[eqleft+varw xbot equalw texth],...
      'horizon','center',...
      'string',''' = ',...
      'visible','off');
   
   ud.h.xder=uicontrol('pos',[eqleft+varw+equalw xbot eqlength texth],...
      'string',ud.o.xder,...
      'horizon','left','style','edit',...
      'call',xder,'visible','off');
   
   ud.h.yvar=uicontrol('pos',[eqleft ybot varw texth],...
      'style','edit',...
      'horizon','right',...
      'string',ud.o.yvar,...
      'call',yname,'visible','off');
   
   eq(3) = uicontrol('style','text',...
      'pos',[eqleft+varw ybot equalw texth],...
      'horizon','center','string',''' = ','visible','off');
   
   ud.h.yder=uicontrol('pos',[eqleft+varw + equalw ybot eqlength texth],...
      'string',ud.o.yder,...
      'horizon','left','style','edit',...
      'call',yder,'visible','off');
   
   
   frame(2) = uicontrol('style','frame','pos',disfrwind,'visible','off');
   
   w1 = [
      'ud = get(gcf,''user'');'...
         'nnn = 	str2num(get(ud.h.wind(1),''string''));'...
         'if isempty(nnn),',...
         '  set(ud.h.wind(1),''string'',''?'');',...
         '  nnn = NaN;',...
         'end,',...
         'ud.c.wind(1) = nnn;',...	    
         'ud.c.name = '''';',...
         'set(gcf,''user'',ud);'];	
   
   w2 = [
      'ud = get(gcf,''user'');'...
         'nnn = 	str2num(get(ud.h.wind(2),''string''));'...
         'if isempty(nnn),',...
         '  set(ud.h.wind(2),''string'',''?'');',...
         '  nnn = NaN;',...
         'end,',...
         'ud.c.wind(2) = nnn;',...	    
         'ud.c.name = '''';',...
         'set(gcf,''user'',ud);'];	
   
   w3 = [
      'ud = get(gcf,''user'');'...
         'nnn = 	str2num(get(ud.h.wind(3),''string''));'...
         'if isempty(nnn),',...
         '  set(ud.h.wind(3),''string'',''?'');',...
         '  nnn = NaN;',...
         'end,',...
         'ud.c.wind(3) = nnn;',...	    
         'ud.c.name = '''';',...
         'set(gcf,''user'',ud);'];	
   
   w4 = [
      'ud = get(gcf,''user'');'...
         'nnn = 	str2num(get(ud.h.wind(4),''string''));'...
         'if isempty(nnn),',...
         '  set(ud.h.wind(4),''string'',''?'');',...
         '  nnn = NaN;',...
         'end,',...
         'ud.c.wind(4) = nnn;',...	    
         'ud.c.name = '''';',...
         'set(gcf,''user'',ud);'];	
   
   
   winbot1 = disfrbot + disfrht - 5 - separation;
   winbot2 = winbot1 - separation;
   winbot3 = winbot2 - separation;
   winbot4 = winbot3 - separation;
   winbot5 = winbot4 - separation;
   
   dwind = uicontrol('style','text',...
      'pos',[eqleft winbot1 disfrw-10 texth],...
      'horizon','center',...
      'string','The display window.','visible','off');
   
   % ud.h.twind contains the handles to the text windows, and ud.h.wind
   % contains the handles to the edit windows.
   
   ud.h.twind(1) = uicontrol('style','text',...
      'pos',[eqleft winbot2 winstrlen texth],...
      'horizon','right',...
      'string',['The minimum value of ',ud.o.xvar,' = '],'visible','off');
   
   ud.h.wind(1) = uicontrol('style','edit',...
      'pos',[eqleft+winstrlen winbot2 40*rr texth],...
      'string',num2str(ud.o.wind(1)),...
      'call',w1,'visible','off');
   
   ud.h.twind(2) = uicontrol('style','text',...
      'pos',[eqleft winbot3 winstrlen texth],...
      'horizon','right',...
      'string',['The maximum value of ',ud.o.xvar,' = '],'visible','off');
   
   ud.h.wind(2) = uicontrol('style','edit',...
      'pos',[eqleft+winstrlen winbot3 40*rr texth],...
      'string',num2str(ud.o.wind(2)),...
      'call',w2,'visible','off');
   
   ud.h.twind(3)= uicontrol('style','text',...
      'pos',[eqleft winbot4 winstrlen texth],...
      'horizon','right',...
      'string',['The minimum value of ',ud.o.yvar,' = '],'visible','off');
   
   ud.h.wind(3) = uicontrol('style','edit',...
      'pos',[eqleft+winstrlen winbot4 40*rr texth],...
      'string',num2str(ud.o.wind(3)),...
      'call',w3,'visible','off');
   
   ud.h.twind(4) = uicontrol('style','text',...
      'pos',[eqleft winbot5 winstrlen texth],...
      'horizon','right',...
      'string',['The maximum value of ',ud.o.yvar,' = '],'visible','off');
   
   ud.h.wind(4) = uicontrol('style','edit',...
      'pos',[eqleft+winstrlen winbot5 40*rr texth],...
      'string',num2str(ud.o.wind(4)),...
      'call',w4,'visible','off');
   
   
   frame(3)=uicontrol('style','frame','pos',pfrwind,'visible','off');
   
   p1 = [
      'ud = get(gcf,''user'');'...
         'ud.c.pname1=get(ud.h.pname1,''string'');'...
         'ud.c.pval1=get(ud.h.pval1,''string'');'...
         'ud.flag = 0;'...
         'ud.c.name = '''';',...
         'set(gcf,''user'',ud);'];
   
   p2 = [
      'ud = get(gcf,''user'');'...
         'ud.c.pname2=get(ud.h.pname2,''string'');'...
         'ud.c.pval2=get(ud.h.pval2,''string'');'...
         'ud.flag = 0;'...
         'ud.c.name = '''';',...
         'set(gcf,''user'',ud);'];
   
   p3 = [
      'ud = get(gcf,''user'');'...
         'ud.c.pname3=get(ud.h.pname3,''string'');'...
         'ud.c.pval3=get(ud.h.pval3,''string'');'...
         'ud.flag = 0;'...
         'ud.c.name = '''';',...
         'set(gcf,''user'',ud);'];
   
   p4 = [
      'ud = get(gcf,''user'');'...
         'ud.c.pname4=get(ud.h.pname4,''string'');'...
         'ud.c.pval4=get(ud.h.pval4,''string'');'...
         'ud.flag = 0;'...
         'ud.c.name = '''';',...
         'set(gcf,''user'',ud);'];
   
   
   pnamew = 40*rr;
   pvalw = 50*rr;
   peqw = 10*rr;
   pbot2 = pfrbot + 5;
   pbot1 = pbot2 + separation;
   pleft1 = eqleft + 50*rr + 5;
   peqleft1 = pleft1 + pnamew;
   pvleft1 = peqleft1 + peqw;
   pleft2 = dfigwidth - 10 - pnamew - pvalw - peqw;
   peqleft2 = pleft2 + pnamew;
   pvleft2 = peqleft2 + peqw;	
   
   paratit=uicontrol('style','text',...
      'pos',[eqleft pfrbot+2 50*rr 2.1*texth],...
      'horizon','center',...
      'string',{'Parameters';'or';'expressions'},'visible','off');
   
   ud.h.pname1 = uicontrol('style','edit',...
      'pos',[pleft1 pbot1 pnamew texth],...
      'horizon','right','string',ud.o.pname1,...
      'call',p1,'visible','off');
   
   equal(1) = uicontrol('style','text',...
      'pos',[peqleft1 pbot1 peqw texth],...
      'horizon','center','string','=','visible','off');
   
   ud.h.pval1 = uicontrol('style','edit',...
      'pos',[pvleft1 pbot1 pvalw texth],...
      'string',ud.o.pval1,'call',p1,'visible','off');
   
   ud.h.pname2 = uicontrol('style','edit',...
      'pos',[pleft2 pbot1 pnamew texth],...
      'horizon','right','string',ud.o.pname2,...
      'call',p2,'visible','off');
   
   equal(2) = uicontrol('style','text',...
      'pos',[peqleft2 pbot1 peqw texth],...
      'horizon','center','string','=','visible','off');
   
   ud.h.pval2 = uicontrol('style','edit',...
      'pos',[pvleft2 pbot1 pvalw texth],...
      'string',ud.o.pval2,'call',p2,'visible','off');
   
   ud.h.pname3 = uicontrol('style','edit',...
      'pos',[pleft1 pbot2 pnamew texth],...
      'horizon','right','string',ud.o.pname3,...
      'call',p3,'visible','off');
   
   equal(3) = uicontrol('style','text',...
      'pos',[peqleft1 pbot2 peqw texth],...
      'horizon','center','string','=','visible','off');
   
   ud.h.pval3 = uicontrol('style','edit',...
      'pos',[pvleft1 pbot2 pvalw texth],...
      'string',ud.o.pval3,...
      'call',p3,'visible','off');
   
   
   ud.h.pname4 = uicontrol('style','edit',...
      'pos',[pleft2 pbot2 pnamew texth],...
      'horizon','right','string',ud.o.pname4,...
      'call',p4,'visible','off');
   
   equal(4) = uicontrol('style','text',...
      'pos',[peqleft2 pbot2 peqw texth],...
      'horizon','center','string','=','visible','off');
   
   ud.h.pval4 = uicontrol('style','edit',...
      'pos',[pvleft2 pbot2 pvalw texth],...
      'string',ud.o.pval4,'call',p4,'visible','off');
   
   butt(1) = uicontrol('style','push',...
      'pos',qwind,...
      'string','Quit','call','pplane5(''quit'')','visible','off');
   
   butt(2) = uicontrol('style','push',...
      'pos',rwind,...
      'string','Revert','call','pplane5(''revert'')','visible','off');
   
   butt(3) = uicontrol('style','push',...
      'pos',pwind,...
      'string','Proceed','call','pplane5(''proceed'')','visible','off');
   
   fframe = uicontrol('style','frame','pos',ffrwind,'visible','off');
   
   
   
   ffrtitle = uicontrol('style','text','pos',[ffrleft+5,winbot1,ffrw-10,texth],...
      'string','The direction field.','horizon','center','visible','off');
   
   radleft = ffrleft + 3;
   radw = 50*rr;
   
   typewindw = radw +6;
   typewind = [ffrleft, ffrbot, typewindw, ffrht-separation-3];
   textwindl = ffrleft+typewindw;
   textleft = textwindl + 3;
   textw = ffrw - typewindw;
   textwind = [textwindl, ffrbot,textw, ffrht-separation-3];
   typeframe = uicontrol('style','frame','pos',typewind,'visible','off');
   
   textframe = uicontrol('style','frame','pos',textwind,'visible','off');
   
   switch ud.o.fieldtype
   case 'nullclines'
      rval1 = 1;rval2 = 0;rval3 = 0;rval4 = 0;
   case 'lines'
      rval1 = 0;rval2 = 2;rval3 = 0;rval4 = 0;
   case 'arrows' 
      rval1 = 0;rval2 = 0;rval3 = 3;rval4 = 0;
   case 'none' 
      rval1 = 0;rval2 = 0;rval3 = 0;rval4 = 4;
   otherwise
      error(['Unknown fieldtype ',ud.o.fieldtype,'.'])
   end		
   
   ud.h.rad(1) = uicontrol('style','radio',...
      'pos',[radleft winbot4 radw texth],...
      'string','Nullclines','value',rval1,'visible','off');
   
   ud.h.rad(2) = uicontrol('style','radio',...
      'pos',[radleft winbot3 radw texth],...
      'string','Lines','value',rval2,'max',2,'visible','off');
   
   ud.h.rad(3) = uicontrol('style','radio',...
      'pos',[radleft winbot2 radw texth],...
      'string','Arrows','value',rval3,'max',3,'visible','off');
   
   ud.h.rad(4) = uicontrol('style','radio',...
      'pos',[radleft winbot5 radw texth],...
      'string','None','value',rval4,'max',4,'visible','off');
   
   for i=1:4
      set(ud.h.rad(i),'user',ud.h.rad(:,[1:(i-1),(i+1):4]));
   end
   
   callrad = [
      'me = get(gcf,''currentobject'');',...
         'kk = get(me,''max'');',...
         'set(get(me,''user''),''value'',0),',...
         'set(me,''value'',kk);',...
         'ud = get(gcf,''user'');',...
         'switch kk,',...
         '       case 1, ud.c.fieldtype = ''nullclines'';',...
         '       case 2, ud.c.fieldtype = ''lines'';',...
         '       case 3, ud.c.fieldtype = ''arrows'';',...
         '       case 4, ud.c.fieldtype = ''none'';',...
         'end,',...
         'set(gcf,''user'',ud);'];
   
   set(ud.h.rad,'call',callrad);
   
   nfptsstr = {'Number of'; 'field points per'; 'row or column.'};
   nfptstext = uicontrol('style','text',...
      'pos',[textleft winbot4 textw-5 2.5*texth],...
      'string',nfptsstr,'horizon','center','visible','off');
   
   callnfpts = [
      'ppset = findobj(''name'',''PPLANE5 Setup'');',...
         'ud = get(ppset,''user'');'...
         'me = ud.h.npts;',...
         'kk = str2num(get(me,''string''));',...
         'if isempty(kk),',...
         '  set(me,''string'',''?'');',...
         '  kk = NaN;',...
         'else,',...
         '  kk = floor(kk);',...
         '  [m,N] = computer;'...
         '  if (N <= 8192),',...
         '	N = 32;',...
         '  else,',...
         '	N = 50;',...
         '  end,'...
         '  kk = min([N,max([5,kk])]);'...
         '  set(me,''string'',num2str(kk));'...
         'end,'...
         'ud.c.npts = kk;',...
         'set(ppset,''user'',ud);'];
   
   ud.h.npts = uicontrol('style','edit',...
      'pos',[textleft+(textw -30*rr)/2,winbot5 30*rr,texth],...
      'string',ud.o.npts,'call',callnfpts,'visible','off');
   
   delgall = ['sud = get(gcf,''user'');',...
         'mh = get(sud.h.gallery,''children'');',...
         'add = findobj(sud.h.gallery,''tag'',''add system'');',...
         'mh(find(mh == add)) = [];',...
         'delete(mh);',...
         'set(sud.h.gallery,''user'',[]);',...
         'set(findobj(''tag'',''load default''),''enable'',''on'')'];
   
   mefile = findobj(ppset,'label','File');
   
   meload = uimenu(mefile,'label','Load a system ...',...
      'call','pplane5(''loadsyst'',''system'');');
   mesave = uimenu(mefile,'label','Save the current system ...',...
      'call','pplane5(''savesyst'',''system'');');
   meloadg = uimenu(mefile,'label','Load a gallery ...',...
      'call','pplane5(''loadsyst'',''gallery'');','separator','on');
   mesaveg = uimenu(mefile,'label','Save a gallery ...',...
      'call','pplane5(''savesyst'',''gallery'');',...
      'tag','savegal');
   medelg = uimenu(mefile,'label','Delete the current gallery',...
      'call',delgall);
   melddg = uimenu(mefile,'label','Load the default gallery',...
      'call','pplane5(''loadsyst'',''default'');','enable','off',...
      'tag','load default');
   meedit = findobj(gcf,'label','Edit'); % uimenu('label','Edit');
   meclrf = uimenu(meedit,'label','Clear equations',...
      'call',['ud = get(gcf,''user'');h = ud.h;',...
         'set([h.xvar,h.xder,h.yvar,h.yder],''string'','''');'],...
      'accelerator','E');
   meclrp = uimenu(meedit,'label','Clear parameters',...
      'call',['ud = get(gcf,''user'');h = ud.h;',...
         'set([h.pname1,h.pname2,h.pname3,h.pname4,h.pval1,h.pval2,h.pval3,h.pval4],''string'','''');'],...
      'accelerator','N');
   meclrwind = uimenu(meedit,'label','Clear display window',...
      'call',['ud = get(gcf,''user'');',...
         'set(ud.h.wind,''string'','''');'],...
      'accelerator','D');
   meclrall = uimenu(meedit,'label','Clear all',...
      'call',['ud = get(gcf,''user'');h = ud.h;',...
         'set(h.wind,''string'','''');',...
         'set([h.pname1,h.pname2,h.pname3,h.pname4,h.pval1,h.pval2,h.pval3,h.pval4],''string'','''');',...
         'set([h.xvar,h.xder,h.yvar,h.yder],''string'','''');'],...
      'accelerator','A',...
      'separator','on');
   merestart = uimenu(mefile,'label','Restart PPLANE5','call',...
      'pplane5(''restart'')','separator','on');
   merevert = uimenu(mefile,'label','Revert','call',...
      'pplane5(''revert'')','separator','on');
   meproceed = uimenu(mefile,'label','Proceed','call',...
      'pplane5(''proceed'')','separator','on',...
      'accelerator','G');
   mequit = uimenu(mefile,'label','Quit PPLANE5','call',...
      'pplane5(''quit'')','separator','on');
   sysmenu = uimenu('label','Gallery','visible','off');
   
   meadd = uimenu(sysmenu,'label','Add current system to the gallery',...
      'call','pplane5(''addgall'');','tag','add system');
   sep = 'on';
   for kk = 1:length(system)
      kkk = num2str(kk);
      if kk == 2, sep = 'off';end
      sysmen(kk) = uimenu(sysmenu,'label',system(kk).name,...
         'call',['pplane5(''system'',',kkk,')'],...
         'separator',sep,'visible','off');
   end
   set(sysmenu,'user',system);
   ud.h.gallery = sysmenu;
   ud.flag = 0;
   ud.egg = (exist('EASTEREGG') ==2);
   
   % Record the handles in the User Data of the Set Up figure.
   
   set(ppset,'user',ud);
   hhhh = findobj(ppset,'type','uicontrol');
   set(hhhh,'units','normal')
   
   set(ppset,'visible','on','resize','on');
   set(get(ppset,'children'),'visible','on');
   set(get(sysmenu,'children'),'visible','on');
 %end
  
elseif strcmp(action,'savesyst')
   
   ppset = findobj('name','PPLANE5 Setup');
   type = input1;
   sud = get(ppset,'user');
   
   switch type
   case 'system'  
      systems = get(sud.h.gallery,'user');
      
      newsyst = sud.c;
      fn = newsyst.name;
      if ~isempty(fn)
         fn(find(abs(fn)==32))='_';   % Replace spaces by underlines.
      end	
      comp = computer;
      switch  comp
      case 'PCWIN'
         filter = [sud.ppdir,'\*.pps'];
      case 'MAC2'
         filter = [sud.ppdir,':*.pps'];
      otherwise
         filter = [sud.ppdir,'/*.pps'];
      end
      [fname,pname] = uiputfile(filter,'Save the system as:');
      if fname == 0,return;end
      if ~strcmp(fname,fn)
         ll = length(fname);
         if (ll>4 & strcmp(fname(ll-3:ll),'.pps'))
            fn = fname(1:(ll-4));
         else
            fn =fname;
         end
         newsyst.name = fn;
         sud.c.name = fn;
         set(ppset,'user',sud);
      end
      newsysts = newsyst;
   case 'gallery'
      systems = get(sud.h.gallery,'user');
      ll = length(systems); 
      if ll == 0
         warndlg(['There are no systems to make up a gallery.'],'Warning');
         return
      end
      names = cell(ll,1);
      for j=1:ll
         names{j} = systems(j).name;
      end	
      [sel,ok] = listdlg('PromptString','Select the systems',...
         'Name','Gallery selection',...
         'ListString',names);
      if isempty(sel)
         return
      else
         newsysts = systems(sel);
      end	
      comp = computer;
      switch  comp
      case 'PCWIN'
         prompt = [sud.ppdir,'\*.ppg'];
      case 'MAC2'
         prompt = [sud.ppdir,':*.ppg'];
      otherwise
         prompt = [sud.ppdir,':/.ppg'];
      end
      [fname,pname] = uiputfile(prompt,'Save the gallery as:');
      
   end  % switch type
   
   ll = length(newsysts);
   newsystsc = cell(19*ll,1);  
   for j = 1:ll
      newsyst = newsysts(j);
      newsyst = struct2cell(newsyst);
      newsyst{15} = num2str(newsyst{15});
      wind = newsyst{16};
      for k = 1:4
         newsyst{15+k} = num2str(wind(k));
      end
      newsystc([(19*(j-1)+1):19*j]) = newsyst;
   end  
   fid = fopen([pname fname],'w');
   for kk = 1:length(newsystc)
      str = newsystc{kk};
      l = length(str);
      if l
         for (k=fliplr(find(str=='\')))
            str = [str(1:k-1) '\' str(k:l)];
            l = l+1;
         end
      end
      fprintf(fid,[str,'\n']);
   end
   fclose(fid);
   
elseif strcmp(action,'loadsyst')  % This loads either a system or a gallery.
   
   sud = get(gcf,'user');
   pos = get(gcf,'pos');
   wpos = [pos(1),pos(2)+pos(4)+20,300,20];
   waith = figure('pos',wpos,...
      'numb','off',...
      'vis','off',...
      'next','replace',...
      'menubar','none',...
      'resize','off',...
      'createfcn','');
   axes('pos',[0.01,0.01,0.98,0.98],...
      'vis','off');
   xp = [0 0 0 0];
   yp = [0 0 1 1];
   xl = [1 0 0 1 1];
   yl = [0 0 1 1 0];
   patchh = patch(xp,yp,'r','edgecolor','r','erase','none');
   lineh = line(xl,yl,'erase','none','color','k');
   type = input1;  
   set(sud.h.gallery,'enable','off');
   if strcmp(type,'default')
      set(waith,'name','Loading the default gallery.','vis','on');
      set(findobj('tag','load default'),'enable','off');
      megall = sud.h.gallery;
      mh = get(megall,'children');
      add = findobj(megall,'tag','add system');
      mh(find(mh == add)) = [];
      delete(mh);
      newsysstruct = get(megall,'user');
      system = sud.system;
      ll = length(system);
      x = 1/(ll+2);
      xp = [xp(2),x,x,xp(2)];
      set(patchh,'xdata',xp);
      set(lineh,'xdata',xl);
      drawnow;
      sep = 'on';
      for kk = 1:length(system)
         kkk = num2str(kk);
         if kk ==2, sep = 'off';end
         uimenu(megall,'label',system(kk).name,...
            'call',['pplane5(''system'',',kkk,')'],...
            'separator',sep);
      end % for
      set(megall,'user',system);
   else
      comp = computer;
      switch  comp
      case 'PCWIN'
         prompt = [sud.ppdir,'\'];
      case 'MAC2'
         prompt = [sud.ppdir,':'];
      otherwise
         prompt = [sud.ppdir,'/'];
      end
      if strcmp(type,'system')
         prompt = [prompt,'*.pps'];
         
         [fname,pname] = uigetfile(prompt,'Select a system to load.');
      elseif strcmp(type,'gallery')
         prompt = [prompt,'*.ppg'];
         
         [fname,pname] = uigetfile(prompt,'Select a gallery to load.');
      end  % if strcmp
      if fname == 0
         delete(waith);
         set(sud.h.gallery,'enable','on');
         return;
      end
      set(waith,'name',['Loading ',fname],'vis','on');
      fid = fopen([pname fname],'r');
      newsysts = {};
      kk = 0;
      while ~feof(fid)
         kk = kk + 1;
         newsysts{kk} = fgetl(fid);
      end
      fclose(fid);
      x = 19/(kk+38);
      xp = [xp(2),x,x,xp(2)];
      set(patchh,'xdata',xp);
      set(lineh,'xdata',xl);
      drawnow;
      newsysts = newsysts([1:kk]);
      if mod(kk,19)
         if strcmp(type,'system')
            warndlg(['The file ',fname, ' does not define a proper system.'],...
               'Warning');
         elseif strcmp(type,'gallery')
            warndlg(['The file ',fname, ' does not define a proper gallery.'],...
               'Warning');
         end
         set(sud.h.gallery,'enable','on');
         delete(waith);
         return
      end % if mod(kk,19)
      nnn = kk/19;
      flds = fields(sud.c);flds=flds(:);
      for j = 0:(nnn-1)
         newsyst = newsysts([(j*19+1):(j+1)*19]);
         wind = newsyst(16:19);
         wind = [str2num(wind{1}),str2num(wind{2}),...
               str2num(wind{3}),str2num(wind{4})];
         newsyst{15} = str2num(newsyst{15});
         newsyst{16} = wind;
         newsyst = newsyst(1:16);
         newsyst = newsyst(:);
         newsysstruct(j+1) = cell2struct(newsyst,flds);
      end  % for j = 0:(nnn-1)	  
   end  % if strcmp(type,'default') & else
   nnn = length(newsysstruct);
   ignoresyst = {};
   for j = 1:nnn
      x = (j+1)/(nnn+2);
      xp = [xp(2),x,x,xp(2)];
      set(patchh,'xdata',xp);
      set(lineh,'xdata',xl);
      drawnow;
      newsyst = newsysstruct(j);
      sname = newsyst.name;
      sname(find(abs(sname) == 95)) = ' '; % Replace underscores with spaces.
      newsyst.name = sname;
      ignore = pplane5('addgall',newsyst);
      if ignore == -1;
         ignoresyst{length(ignoresyst)+1} = sname;
      end  
   end % for j = 1:nnn
   l = length(ignoresyst);
   if l  % There was at least one system which was a dup with a
      % different name.
      if l == 1
         message = {['The system "',ignoresyst{1},'" duplicates a ',...
                  'system already in the gallery and was not added.']};
      else
         message = 'The systems ';
         for k = 1:(l-1)
            message = [message,'"',ignoresyst{k},'", ']; 
         end
         message = {[message,'and "',ignoresyst{l},'" duplicate ',...
                  'systems already in the gallery and were not added.']};
      end % if l == 1 & else
      helpdlg(message,'Ignored systems');
   end  % if l
   if strcmp(type,'system') % Added a system.
      if ignore > 0 % The system was ignored.
         kk = ignore;
      else
         systems = get(sud.h.gallery,'user');
         kk = length(systems);
      end
      pplane5('system',kk);
   end  
   if strcmp('type','default')
      pplane5('system',1);
   end  
   set(sud.h.gallery,'enable','on');
   x = 1;
   xp = [xp(2),x,x,xp(2)];
   set(patchh,'xdata',xp);
   set(lineh,'xdata',xl);
   drawnow;
   delete(waith);
   
elseif strcmp(action,'addgall')
   
   output = 0;
   ppset = findobj('name','PPLANE5 Setup');
   sud = get(ppset,'user');
   if nargin < 2    % We are adding the current system.
      
      syst = sud.c;
      sname = inputdlg('Provide a name for this system.','System name',1,{syst.name});
      if isempty(sname),return;end
      sname = sname{1};
      if ~strcmp(sname,syst.name)
         sud.c.name = sname;
         set(ppset,'user',sud);
         syst.name = sname;
      end  
   else  % We have a system coming from a file.
      syst = input1;
      sname = syst.name;
   end
   systems = get(sud.h.gallery,'user');
   ll = length(systems);
   kk = 1;
   while ((kk<=ll) & (~strcmp(sname,systems(kk).name)))
      kk = kk + 1;
   end
   nameflag = (kk<=ll);
   ssyst = rmfield(syst,'name');
   kk = 1;
   while ((kk<=ll) & (~isequal(ssyst,rmfield(systems(kk),'name'))))
      kk = kk + 1;
   end
   systflag = 2*(kk<=ll);
   flag = nameflag + systflag;
   switch flag
   case 1  % Same name but different system.
      mh = findobj(sud.h.gallery,'label',sname);
      prompt = {['The system "',sname,'", which you wish to add to the gallery has ',...
               'the same name as a different system already in the gallery.  Please ',...
               'specify the name for the newly added system.'],...
            'Specify the name for the old system.'};
      title = 'Two systems with the same name';
      lineno = 1;
      defans = {sname,sname};
      answer = inputdlg(prompt,title,lineno,defans);
      if isempty(answer),break,end
      sname = answer{1};
      systems(kk).name = answer{2};
      set(mh,'label',answer{2});
      output = kk;
   case 2  % Two names for the same system.
      oldname = systems(kk).name;
      mh = findobj(sud.h.gallery,'label',oldname);
      
      prompt = {['The system "',sname,'", which you wish to add ',...
               'to the gallery is the same as a system which is ',...
               'already in the gallery with the name "',oldname,'".  ',...
               'Please specify which name you wish to use.']};
      title = 'Two names for the same system.';
      lineno = 1;
      defans = {oldname};
      answer = inputdlg(prompt,title,lineno,defans);
      if isempty(answer),break,end
      systems(kk).name = answer{1};
      set(mh,'label',answer{1});
      output = kk;
   case 3 % Systems and names the same.
      output = -1;
   otherwise
   end  % switch
   set(sud.h.gallery,'user',systems);
   syst.name = sname;
   if flag <=1
      switch ll
      case 0
         systems = syst;
         sepstr = 'on';
      case 8
         systems(9) = syst;
         if strcmp(systems(8).name,'Duffing''s equation')
            sepstr = 'on';
         else
            sepstr = 'off';
         end
      otherwise
         systems(ll+1) = syst;
         sepstr = 'off';
      end      
      kkk = num2str(ll+1); 
      newmenu = uimenu(sud.h.gallery,'label',sname,...
         'call',['pplane5(''system'',',kkk,')'],...
         'separator',sepstr);
      set(findobj('tag','savegal'),'enable','on');
   end
   set(sud.h.gallery,'user',systems);
   
elseif strcmp(action,'system')
   
   ppset = findobj('name','PPLANE5 Setup');
   ud = get(ppset,'user');
   kk = input1;
   if isstr(kk)
      kk = str2num(input1);
   end
   system = get(ud.h.gallery,'user');
   syst = system(kk);
   ud.o = syst;
   ud.c = syst;
   xname = syst.xvar;
   yname = syst.yvar;
   set(ud.h.xvar,'string',xname);
   set(ud.h.yvar,'string',yname);
   set(ud.h.xder,'string',syst.xder);
   set(ud.h.yder,'string',syst.yder);
   set(ud.h.pname1,'string',syst.pname1);
   set(ud.h.pname2,'string',syst.pname2);
   set(ud.h.pname3,'string',syst.pname3);
   set(ud.h.pname4,'string',syst.pname4);
   set(ud.h.pval1,'string',syst.pval1);
   set(ud.h.pval2,'string',syst.pval2);
   set(ud.h.pval3,'string',syst.pval3);
   set(ud.h.pval4,'string',syst.pval4);
   set(ud.h.twind(1),'string',['The minimum value of ',xname,' = ']);
   set(ud.h.twind(2),'string',['The maximum value of ',xname,' = ']);
   set(ud.h.twind(3),'string',['The minimum value of ',yname,' = ']);
   set(ud.h.twind(4),'string',['The maximum value of ',yname,' = ']);
   set(ud.h.wind(1),'string',num2str(syst.wind(1)));
   set(ud.h.wind(2),'string',num2str(syst.wind(2)));
   set(ud.h.wind(3),'string',num2str(syst.wind(3)));
   set(ud.h.wind(4),'string',num2str(syst.wind(4)));
   set(ud.h.npts,'string',num2str(syst.npts));
   
   switch syst.fieldtype
   case 'nullclines'
      rval = [1 0 0 0];
   case 'lines'
      rval = [0 2 0 0];
   case 'arrows' 
      rval = [0 0 3 0];
   case 'none' 
      rval = [0 0 0 4];
   otherwise
      error(['Unknown fieldtype ',ud.o.fieldtype,'.'])
   end		
   
   for i=1:4
      set(ud.h.rad(i),'value',rval(i));
   end	
   ud.flag = 0;
   set(ppset,'user',ud);
   
elseif strcmp(action,'revert')
   
   ud = get(gcf,'user');
   ud.c = ud.o;
   syst = ud.o;
   xname = syst.xvar;
   yname = syst.yvar;
   set(ud.h.xvar,'string',xname);
   set(ud.h.yvar,'string',yname);
   set(ud.h.xder,'string',syst.xder);
   set(ud.h.yder,'string',syst.yder);
   set(ud.h.pname1,'string',syst.pname1);
   set(ud.h.pname2,'string',syst.pname2);
   set(ud.h.pname3,'string',syst.pname3);
   set(ud.h.pname4,'string',syst.pname4);
   set(ud.h.pval1,'string',syst.pval1);
   set(ud.h.pval2,'string',syst.pval2);
   set(ud.h.pval3,'string',syst.pval3);
   set(ud.h.pval4,'string',syst.pval4);
   set(ud.h.twind(1),'string',['The minimum value of ',xname,' = ']);
   set(ud.h.twind(2),'string',['The maximum value of ',xname,' = ']);
   set(ud.h.twind(3),'string',['The minimum value of ',yname,' = ']);
   set(ud.h.twind(4),'string',['The maximum value of ',yname,' = ']);
   set(ud.h.wind(1),'string',num2str(syst.wind(1)));
   set(ud.h.wind(2),'string',num2str(syst.wind(2)));
   set(ud.h.wind(3),'string',num2str(syst.wind(3)));
   set(ud.h.wind(4),'string',num2str(syst.wind(4)));
   set(ud.h.npts,'string',num2str(syst.npts));
   
   switch syst.fieldtype
   case 'lines'
      rval(1) = 1;
      rval(2) = 0;rval(3) = 0;
   case 'arrows' 
      rval(1) = 0;rval(2) = 2;rval(3) = 0;
   case 'none' 
      rval(1) = 0;rval(2) = 0;rval(3) = 3;
   otherwise
      error(['Unknown fieldtype ',ud.o.fieldtype,'.'])
   end		
   
   for i=1:3
      set(ud.h.rad(i),'value',rval(i));
   end	
   set(gcf,'user',ud);
   
elseif strcmp(action,'proceed')
   
   % Proceed connects Setup with the Display window.
   
   ppset = gcf;
   sud = get(ppset,'user');
   sud.o = sud.c;
   set(ppset,'user',sud);
   he = findobj('name','PPLANE5 Equilibrium point data');
   hl = findobj('name','PPLANE5 Linearization');
   close([he;hl]);
   % set([he;hl],'vis','off');
   
   % Some error checking that has to be done no matter what.
   
   WINvect = sud.c.wind;
   if any(isnan(WINvect))
      sud.flag = 0;
      set(ppset,'user',sud);
      errmsg = ['One of the entries defining the display window ',...
            'is not a number.'];
      fprintf('\a')
      errordlg(errmsg,'PPLANE error','on');
      return
   end
   xstr = sud.c.xvar;
   if isempty(xstr)
      sud.flag = 0;
      set(ppset,'user',sud);
      errmsg = 'The first dependent variable needs a name.';
      fprintf('\a')
      errordlg(errmsg,'PPLANE error','on');
      return
   end
   ystr = sud.c.yvar;
   if isempty(ystr)
      sud.flag = 0;
      set(ppset,'user',sud);
      errmsg = 'The second dependent variable needs a name.';
      fprintf('\a')
      errordlg(errmsg,'PPLANE error','on');
      return
   end
   if WINvect(2)<= WINvect(1)
      sud.flag = 0;
      set(ppset,'user',sud);
      errmsg = ['The minimum value of ', xstr,...
            ' must be smaller than the maximum value.'];
      fprintf('\a')
      errordlg(errmsg,'PPLANE error','on');
      return
   end
   if WINvect(4)<= WINvect(3)
      sud.flag = 0;
      set(ppset,'user',sud);
      errmsg = ['The minimum value of ', ystr,...
            ' must be smaller than the maximum value.'];
      fprintf('\a')
      errordlg(errmsg,'PPLANE error','on');
      return
   end
   if isnan(sud.c.npts)
      sud.flag = 0;
      set(ppset,'user',sud);
      errmsg = 'The entry for the number of field points is not a number.';
      fprintf('\a')
      errordlg(errmsg,'PPLANE error','on');
      return
   end
   
   %  sud.flag = 0 if this is the first time through for this equation, 
   %  sud.flag = 1 if only the window dimensions or the field data 
   %  have been changed.
   
   % If sud.flag == 1 we only have to update things.
   
   if (sud.flag == 1)
      Arrflag = sud.c.fieldtype;
      NumbFPts = sud.c.npts;
      ppdisp = findobj('name','PPLANE5 Display');
      dud = get(ppdisp,'user');
      aud = get(dud.axes,'user');
      wind = sud.c.wind(:);
      if (~all(wind == dud.syst.wind(:)))
         dwind = [wind(1); wind(3); -wind(2); -wind(4)];
         aud.DY = [wind(2)-wind(1); wind(4)-wind(3)];
         aud.cwind = dwind - dud.settings.magn*[aud.DY;aud.DY];
         set(dud.axes,'user',aud);
      end
      
      if ~strcmp(Arrflag,dud.syst.fieldtype)
         arr = dud.arr;
         menull = findobj('tag','null');      
         switch Arrflag
         case 'nullclines'
            set([arr.hx;arr.hy;arr.barrows],'vis','on');
            set([arr.lines;arr.arrows],'vis','off');
            set(menull,'enable','off','label','Plot nullclines.');	  
         case 'lines'
            set(arr.lines,'vis','on');
            set([arr.hx;arr.hy;arr.arrows;arr.barrows],'vis','off');
            set(menull,'enable','on');	  
         case 'arrows'
            set(arr.arrows,'vis','on');
            set([arr.lines;arr.hx;arr.hy;arr.barrows],'vis','off');
            set(menull,'enable','on');	  
         otherwise
            set([arr.hx;arr.hy;arr.lines;arr.arrows;arr.barrows],'vis','off');
            set(menull,'enable','on');	  
         end
         dud.syst.fieldtype = Arrflag;
         set(ppdisp,'user',dud);
      end
      
      if ( (NumbFPts ~= dud.syst.npts) | (any(WINvect ~= dud.syst.wind) ) )
         dud.syst.wind = WINvect;			
         dud.syst.npts = NumbFPts;
         set(ppdisp,'user',dud);
         pplane5('dirfield',ppdisp);
      end
   else
      sud.flag = 1;
      set(ppset,'user',sud);	
      
      sud = get(ppset,'user');
      %  WINvect = sud.c.wind;
      
      Xname = sud.c.xvar;
      Yname = sud.c.yvar;
      xderivstr = sud.c.xder;
      yderivstr = sud.c.yder;
      pname1 = sud.c.pname1;
      pname2 = sud.c.pname2;
      pname3 = sud.c.pname3;
      pname4 = sud.c.pname4;
      parav1 = sud.c.pval1;
      parav2 = sud.c.pval2;
      parav3 = sud.c.pval3;
      parav4 = sud.c.pval4;
      
      % Convert the parameters to their current values.
      % First remove the blanks.
      
      xderivstr(find(abs(xderivstr)==32))=[];
      yderivstr(find(abs(yderivstr)==32))=[];
      if ~isempty(parav1)
	parav1(find(abs(parav1)==32))=[];
      end
      if ~isempty(parav2)
	parav2(find(abs(parav2)==32))=[];
      end
      if ~isempty(parav3)
	parav3(find(abs(parav3)==32))=[];
      end
      if ~isempty(parav4)
	parav4(find(abs(parav4)==32))=[];
      end
      
      % Next remove the periods inserted by users attempting to make the 
      % function array smart.
      
      l=length(xderivstr);
      for ( k = fliplr(findstr('.',xderivstr)))
         if (find('*/^' == xderivstr(k+1)))
            xderivstr = [xderivstr(1:k-1), xderivstr(k+1:l)];
         end
         l=l-1;
      end
      
      l=length(yderivstr);
      for ( k = fliplr(findstr('.',yderivstr)))
         if (find('*/^' == yderivstr(k+1)))
            yderivstr = [yderivstr(1:k-1), yderivstr(k+1:l)];
         end
         l=l-1;
       end
       
       l=length(parav1);
       for ( k = fliplr(findstr('.',parav1)))
	 if (find('*/^' == parav1(k+1)))
	   parav1 = [parav1(1:k-1), parav1(k+1:l)];
	 end
	 l=l-1;
       end
    
       l=length(parav2);
       for ( k = fliplr(findstr('.',parav2)))
	 if (find('*/^' == parav2(k+1)))
	   parav2 = [parav2(1:k-1), parav2(k+1:l)];
	 end
	 l=l-1;
       end
        
       l=length(parav3);
       for ( k = fliplr(findstr('.',parav3)))
	 if (find('*/^' == parav3(k+1)))
	   parav3 = [parav3(1:k-1), parav3(k+1:l)];
	 end
	 l=l-1;
       end
       
       l=length(parav4);
       for ( k = fliplr(findstr('.',parav4)))
	 if (find('*/^' == parav4(k+1)))
	   parav4 = [parav4(1:k-1), parav4(k+1:l)];
	 end
	 l=l-1;
       end
        
       
       % Build strings for the title.
      txderstr = xderivstr;
      tyderstr = yderivstr;

      kxder = find(abs(txderstr)==42);
      txderstr(kxder)=' '*ones(size(kxder));
      txderstr = strrep(txderstr,'-',' - ');
      txderstr = strrep(txderstr,'+',' + ');
      if (abs(txderstr(1)) == 32)
	txderstr = txderstr(2:length(txderstr));
      end
       
      kxder = find(abs(tyderstr)==42);
      tyderstr(kxder)=' '*ones(size(kxder));
      tyderstr = strrep(tyderstr,'-',' - ');
      tyderstr = strrep(tyderstr,'+',' + ');
      if (abs(tyderstr(1)) == 32)
	tyderstr = tyderstr(2:length(tyderstr));
      end
      
      tstr1 = [Xname,' '' = ', txderstr];
      tstr2 = [Yname,' '' = ', tyderstr];
      tstr = str2mat(tstr1,tstr2);
      dud.tstr = tstr;
      pstr1 = {' ';' '};
      pstr2 = {' ';' '};
      
      pstring = cell(4,1);
      pstring{1} = [pname1,' = ',parav1];
      pstring{2} = [pname2,' = ',parav2];
      pstring{3} = [pname3,' = ',parav3];
      pstring{4} = [pname4,' = ',parav4];
      
    if ~isempty(parav1)
      tp1str = parav1;
      
      kxder = find(abs(tp1str)==42);  % Get rid of *s
      tp1str(kxder)=' '*ones(size(kxder));
      tp1str = strrep(tp1str,'-',' - ');  % Extra spaces
      tp1str = strrep(tp1str,'+',' + ');
      if (abs(tp1str(1)) == 32)  % Get rid of starting space
	tp1str = tp1str(2:length(tp1str));
      end
      pstring{1} = [pname1,' = ', tp1str];
    end
    
    if ~isempty(parav2)
      tp2str = parav2;
      
      kxder = find(abs(tp2str)==42);  % Get rid of *s
      tp2str(kxder)=' '*ones(size(kxder));
      tp2str = strrep(tp2str,'-',' - ');  % Extra spaces
      tp2str = strrep(tp2str,'+',' + ');
      if (abs(tp2str(1)) == 32)  % Get rid of starting space
	tp2str = tp2str(2:length(tp2str));
      end
      pstring{2} = [pname2,' = ', tp2str];
    end
    
    if ~isempty(parav3)
      tp3str = parav3;
      
      kxder = find(abs(tp3str)==42);  % Get rid of *s
      tp3str(kxder)=' '*ones(size(kxder));
      tp3str = strrep(tp3str,'-',' - ');  % Extra spaces
      tp3str = strrep(tp3str,'+',' + ');
      if (abs(tp3str(1)) == 32)  % Get rid of starting space
	tp3str = tp3str(2:length(tp3str));
      end
      pstring{3} = [pname3,' = ', tp3str];
    end
    
    if ~isempty(parav4)
      tp4str = parav4;
      
      kxder = find(abs(tp4str)==42);  % Get rid of *s
      tp4str(kxder)=' '*ones(size(kxder));
      tp4str = strrep(tp4str,'-',' - ');  % Extra spaces
      tp4str = strrep(tp4str,'+',' + ');
      if (abs(tp4str(1)) == 32)  % Get rid of starting space
	tp4str = tp4str(2:length(tp4str));
      end
      pstring{4} = [pname4,' = ', tp4str];
    end
   
    % Get ready to do some error trapping.
    
    SS = warning;
    warning off
    XxXxXx = WINvect(1) + rand*(WINvect(2)-WINvect(1));
    YyYyYy = WINvect(3) + rand*(WINvect(4)-WINvect(3));
    err = 0;

      % Now we remove the backslashes (\) used to get Greek into the
      % labels. 
      txname = Xname;
      tyname = Yname;    
      xderivstr(find(abs(xderivstr)==92))=[];
      yderivstr(find(abs(yderivstr)==92))=[];
      Xname(find(abs(Xname)==92))=[];
      Yname(find(abs(Yname)==92))=[];
      
      eval([Xname,'=XxXxXx;'],'err = 1;');
      if err
         sud.flag = 0;
         set(ppset,'user',sud);
         errmsg = ['"',xstr, '" is not a valid variable name in MATLAB.'];
         fprintf('\a')
         errordlg(errmsg,'PPLANE error','on');
         return
      end
      err = 0;
      eval([Yname,'=YyYyYy;'],'err = 1;');
      if err
         sud.flag = 0;
         set(ppset,'user',sud);
         errmsg = ['"',ystr, '" is not a valid variable name in MATLAB.'];
         fprintf('\a')
         errordlg(errmsg,'PPLANE error','on');
         return
      end
      
    
    % Replace the parameters/expressions with their values.
    pflag = zeros(1,4);
    perr = [];
    if ~isempty(pname1)
      if isempty(parav1) 
	perr = sud.h.pval1;
      else 
	if isempty(str2num(parav1)) % This is an expression.
	  tparav1 = parav1;
	  tparav1(find(abs(tparav1)==92))=[];
	  err = 0; pval = '';
	  eval(['pval = ',tparav1,';'],'err=1;');
	  if (err | isempty(pval))
	    errmsg = ['The expression for ',pname1,' is not valid.'];
	    fprintf('\a')
	    errordlg(errmsg,'DFIELD5 error','on');
	    warning(SS)
	    return
	  end
	end
	xxderivstr = dfield5('paraeval',pname1,parav1,xderivstr);
	yyderivstr = pplane5('paraeval',pname1,parav1,yderivstr);
	if (~strcmp(xxderivstr,xderivstr) | ~strcmp(yyderivstr,yderivstr) )
	  pflag(1) = 1;
	  xderivstr = xxderivstr;
	  yderivstr = yyderivstr;
	end
      end
    end

    if ~isempty(pname2)
      if isempty(parav2) 
	perr = sud.h.pval2;
      else 
	if isempty(str2num(parav2)) % This is an expression.
	  tparav2 = parav2;
	  tparav2(find(abs(tparav2)==92))=[];
	  err = 0; pval = '';
	  eval(['pval = ',tparav2,';'],'err=1;');
	  if (err | isempty(pval))
	    errmsg = ['The expression for ',pname2,' is not valid.'];
	    fprintf('\a')
	    errordlg(errmsg,'DFIELD5 error','on');
	    warning(SS)
	    return
	  end
	end
	xxderivstr = dfield5('paraeval',pname2,parav2,xderivstr);
	yyderivstr = pplane5('paraeval',pname2,parav2,yderivstr);
	if (~strcmp(xxderivstr,xderivstr) | ~strcmp(yyderivstr,yderivstr) )
	  pflag(2) = 1;
	  xderivstr = xxderivstr;
	  yderivstr = yyderivstr;
	end
      end
    end

    if ~isempty(pname3)
      if isempty(parav3) 
	perr = sud.h.pval3;
      else 
	if isempty(str2num(parav3)) % This is an expression.
	  tparav3 = parav3;
	  tparav3(find(abs(tparav3)==92))=[];
	  err = 0; pval = '';
	  eval(['pval = ',tparav3,';'],'err=1;');
	  if (err | isempty(pval))
	    errmsg = ['The expression for ',pname3,' is not valid.'];
	    fprintf('\a')
	    errordlg(errmsg,'DFIELD5 error','on');
	    warning(SS)
	    return
	  end
	end
	xxderivstr = dfield5('paraeval',pname3,parav3,xderivstr);
	yyderivstr = pplane5('paraeval',pname3,parav3,yderivstr);
	if (~strcmp(xxderivstr,xderivstr) | ~strcmp(yyderivstr,yderivstr) )
	  pflag(3) = 1;
	  xderivstr = xxderivstr;
	  yderivstr = yyderivstr;
	end
      end
    end

    if ~isempty(pname4)
      if isempty(parav4) 
	perr = sud.h.pval4
      else 
	if isempty(str2num(parav4)) % This is an expression.
	  tparav4 = parav4;
	  tparav4(find(abs(tparav4)==92))=[];
	  err = 0; pval = '';
	  eval(['pval = ',tparav4,';'],'err=1;');
	  if (err | isempty(pval))
	    errmsg = ['The expression for ',pname4,' is not valid.'];
	    fprintf('\a')
	    errordlg(errmsg,'DFIELD5 error','on');
	    warning(SS)
	    return
	  end
	end
	xxderivstr = dfield5('paraeval',pname4,parav4,xderivstr);
	yyderivstr = pplane5('paraeval',pname4,parav4,yderivstr);
	if (~strcmp(xxderivstr,xderivstr) | ~strcmp(yyderivstr,yderivstr) )
	  pflag(4) = 1;
	  xderivstr = xxderivstr;
	  yderivstr = yyderivstr;
	end
      end
    end

       
      % We have to make the derivative strings array smart.
      
      l = length(xderivstr);
      for (k=fliplr(find((xderivstr=='^')|(xderivstr=='*')|(xderivstr=='/'))))
         xderivstr = [xderivstr(1:k-1) '.' xderivstr(k:l)];
         l = l+1;
      end
      
      l = length(yderivstr);
      for (k=fliplr(find((yderivstr=='^')|(yderivstr=='*')|(yderivstr=='/'))))
         yderivstr = [yderivstr(1:k-1) '.' yderivstr(k:l)];
         l = l+1;
      end
      
      % Somemore error trapping.
      
      err = 0;res = 1;
      eval(['res = ',xderivstr, ';'],'err = 1;');
      if err | isempty(res)
         if isempty(perr)
            errmsg = ['The first differential equation ',...
                  'is not entered correctly.'];
         else
            errstr1 = ['The first differential equation ',...
                  'does not evaluate correctly.'];
            errstr2 = ['At least one of the parameter values is not ',...
                  'a number.'];
            errmsg = str2mat(errstr1,errstr2);
            set(perr,'string','?');
         end
         sud.flag = 0;
         set(ppset,'user',sud);
         fprintf('\a')
         errordlg(errmsg,'PPLANE error','on');
         return;
       end

      err = 0;res = 1;
      eval(['res = ',yderivstr, ';'],'err = 1;');
      if err | isempty(res)
         if isempty(perr)
            errmsg = ['The second differential equation ',...
                  'is not entered correctly.'];
         else
            errstr1 = ['The second differential equation ',...
                  'does not evaluate correctly.'];
            errstr2 = ['At least one of the parameter values is not ',...
                  'a number.'];
            errmsg = str2mat(errstr1,errstr2);
            set(perr,'string','?');
         end
         sud.flag = 0;
         set(ppset,'user',sud);
         fprintf('\a')
         errordlg(errmsg,'PPLANE error','on');
         return;
      end
      
      % If an old function m-file exists delete it, and then build a new one.
      
      %    if (~strcmp(dfcn,'') & exist(dfcn)==2) delete([dfcn,'.m']);end
      tee = clock;
      tee = ceil(tee(6)*100);
      dfcn=['pptp',num2str(tee)];
      fcnstr = ['function YYyYypr = ',dfcn,'(t,YyYy)\n\n'];
      commstr = '%%%% Created by PPLANE5\n\n';
      varstr = [Xname,' = YyYy(1,:);', Yname,' = YyYy(2,:);\n\n'];
      lenstr = ['l = length(YyYy(1,:));\n'];
      derstr1 = ['XxXxxpr = ', xderivstr,';\n'];
      derstr2 = ['if (length(XxXxxpr) < l) XxXxxpr = XxXxxpr*ones(1,l);end\n'];
      derstr3 = ['YyYyypr = ', yderivstr,';\n'];
      derstr4 = ['if (length(YyYyypr) < l) YyYyypr = YyYyypr*ones(1,l);end\n'];
      derstr5 = 'YYyYypr = [XxXxxpr;YyYyypr];\n';
      ppf = fopen([dfcn,'.m'],'w');
      
      fprintf(ppf,fcnstr);
      fprintf(ppf,commstr);
      fprintf(ppf,varstr);
      fprintf(ppf,lenstr);
      fprintf(ppf,derstr1);
      fprintf(ppf,derstr2);
      fprintf(ppf,derstr3);
      fprintf(ppf,derstr4);
      fprintf(ppf,derstr5);
      fclose(ppf);
      
      
      
      % Find PPLANE5 Display if it exists.	
      % If PPLANE5 Display exists, update it.  If it does not build it.
      
      ppdisp = findobj('name','PPLANE5 Display');
      if (~isempty(ppdisp))
         figure(ppdisp);
         dud = get(ppdisp,'user');
         dud.syst = sud.c;
         dud.settings = sud.settings;
         dfcnn = dud.function;
         if  exist(dfcnn)==2, delete([dfcnn,'.m']);end
         xmstr = [txname,' vs. t'];
         ymstr = [tyname,' vs. t'];
         set(dud.menu(3),'label',xmstr);
         set(dud.menu(5),'label',ymstr);
         menull = findobj('tag','null');
         if ~isempty(menull)
            delete(get(menull,'user'));
            set(menull,'user',[]);      
         end
      else
         ppdisp = figure('name','PPLANE5 Display',...
            'numb','off',...
            'interrupt','on',...
            'visible','off',...
            'tag','PPLANE5');
         pplane5('figdefault',ppdisp);
         dud = get(ppdisp,'user');
         dud.syst = sud.c;
         dud.settings = sud.settings;
         dud.egg = sud.egg;      
         fs = dud.fontsize;
         ssize = dud.ssize;
         r = ssize/10;
         ppaxw = 437*1.2;	% Default axes width
         ppaxh = 315*1.2;	% Default axes height
         ppaxl = 45*1.2;	% Default axes left
         buttw = 40*1.2;     % Default button width
         titleh = 45;    % Default title height.  This is changed later.
         nframeh = 70;   % Default notice frame height
         ppaxb = 4+nframeh+35;
         bottomedge = 38;
         ppdh = bottomedge + nframeh + ppaxh + titleh;
         uni = get(0,'units');
         set(0,'units','pixels');
         ss = get(0,'screensize');
         set(0,'units',uni);
         sw = ss(3);sh = ss(4);
	 bottom = 10;
         if r*ppdh > sh -bottom -35;
            r = (sh-bottom-35)/ppdh;
            fs = 10*r;
            lw = 0.5*r;
            set(gcf,'defaultaxesfontsize',fs,'defaultaxeslinewidth',lw);
            set(gcf,'defaulttextfontsize',fs);
            set(gcf,'defaultlinelinewidth',lw);
            set(gcf,'defaultuicontrolfontsize',fs*0.9);
         end
         
         % Set up the bulletin window.
         
         nframe = uicontrol('style','frame','visible','on');
         nstr = {'More';'than';'five';'lines';'of text'};
         dud.notice = uicontrol('style','text',...
            'horiz','left',...
            'string',nstr,'visible','on');
         ext = get(dud.notice,'extent');
         nframeh = ext(4)+2;
         
         
         titleh = r*titleh;
         ppaxl = r*ppaxl;
         ppaxw = r*ppaxw;
         ppaxh = r*ppaxh;
         ppaxb = nframeh+r*bottomedge;
         buttw = r*buttw;
         butth = fs+10*r;
         buttl = ppaxl + ppaxw + 5;
         buttsep = (ppaxh - butth)/2;
         
         % Set up the plot axes.
         
         dud.axes = axes('units','pix',...
            'position',[ppaxl,ppaxb,ppaxw,ppaxh],...
            'next','add',...
            'box','on',...
            'interrupt','on',...
            'xgrid','on',...
            'ygrid','on',...
            'drawmode','fast',...
            'visible','off',...
            'tag','display axes');
         
         % Set up the title.
         
         dud.title.axes = axes('box','off','xlim',[0 1],'ylim',[0 1],...
            'units','pix','vis','off','xtick',[-1],'ytick',[-1],...
            'xticklabel','','yticklabel','');
         
         dud.title.eq = text(0.07,0.5,' ','vert','middle');
         dud.title.p1 = text(0.75,0.5,' ','vert','middle');
         dud.title.p2 = text(0.65,0.5,' ','vert','middle');
         tstr = {'x_2';'y^2'};
         set(dud.title.eq,'string',tstr,'units','pix');
         ext = get(dud.title.eq ,'extent');
         titleh = ext(4)+15*r;
         set(dud.title.eq,'units','data');
         taxpos = [ppaxl,ppaxb+ppaxh,ppaxw,titleh];
         set(dud.title.axes,'pos',taxpos,'color',get(gcf,'color'));
         
         % Finish the positions.
         
         ppdw = buttl + buttw +5;
         ppdh = ppaxb+ppaxh+titleh;
         set(nframe,'pos',[10,1,ppdw-20,nframeh]);
         set(dud.notice,'pos',[15,1,ppdw-30,nframeh-2],...
            'string',{' ';' ';' ';' ';' '});
         
         ppdleft = max((sw-ppdw)/2,sw-ppdw-60);
         ppdbot = sh - ppdh - 35;
         ppdpos = [ppdleft,ppdbot,ppdw,ppdh];
         set(ppdisp,'resize','on');
         set(ppdisp,'pos',ppdpos);
         
         Arrflag = sud.c.fieldtype;
         
         % Set up the buttons and the menu.
         
         stopstr = 'aud = get(gca,''user'');aud.stop = 4;set(gca,''user'',aud);';
         
         dbutt(1) = uicontrol('style','push',...
            'pos',[buttl,ppaxb+2*buttsep,buttw,butth],...
            'string','Stop','call',stopstr,'vis','off','tag','stop');
         
         dbutt(2) = uicontrol('style','push',...
            'pos',[buttl,ppaxb,buttw,butth],...
            'string','Quit',...
            'call','pplane5(''quit'')','visible','off');
         
         dbutt(3) = uicontrol('style','push',...
            'pos',[buttl,ppaxb+buttsep,buttw,butth],...
            'string','Print',...
            'call','pplane5(''print'')','visible','off');	
         
         dud.butt = dbutt;
         
         fmenu = findobj(gcf,'label','File');  % uimenu('label','File');
         emenu = findobj(gcf,'label','Edit'); % uimenu('label','Edit');
         
         solmenu = uimenu('label','Solutions');
         
         megraph  = uimenu('label','Graph','visible','off');
         menu(1) = uimenu('label','Options','visible','off');
         
         meprint = uimenu(fmenu,'label','Print PPLANE5 Display',...
            'Handlevis','on', ...
            'Accelerator','P', ...
            'Callback','print -noui');
         
         meprsetup =   uimenu(fmenu,'label','Page Position ...',...
            'Handlevis','on', ...
            'Callback','pagedlg(gcf)'); 
         
         meprdlg =   uimenu(fmenu,'label','Print ...',...
            'Handlevis','on', ...
            'Callback','printdlg(gcf)');  
         
         mequit = uimenu(fmenu,'label','Quit PPLANE5',...
            'call','pplane5(''quit'')','separator','on');
         
         menukey = uimenu(solmenu,'label','Keyboard input.','call',...
            'pplane5(''kbd'')');
         
         mesev   = uimenu(solmenu,'label','Plot several solutions.',...
            'call','pplane5(''several'')');
         
         menu(3) = uimenu(megraph,'label',[txname,' vs. t'],...
            'call','pplane5(''plotxy'',1)');
         
         menu(5) = uimenu(megraph,'label',[tyname,' vs. t'],...
            'call','pplane5(''plotxy'',2)');
         
         meplot3 = uimenu(megraph,'label','Both',...
            'call','pplane5(''plotxy'',3)');
         
         meplot4 = uimenu(megraph,'label','3 D',...
            'call','pplane5(''plotxy'',4)');
         
         meplot5 = uimenu(megraph,'label','Composite',...
            'call','pplane5(''plotxy'',5)');
         
         menu(2) = uimenu(emenu,'label','Zoom in.',...
            'call','pplane5(''zoomin'')');
         
         zsqmenu = uimenu(emenu,'label','Zoom in square.',...
            'call','pplane5(''zoominsq'')');
         
         zbmenu = uimenu(emenu,'label','Zoom back.','call',...
            'pplane5(''zoomback'')','enable','off','tag','zbmenu');
         
         
         meeqpt  = uimenu(solmenu,'label','Find an equilibrium point.',...
            'call','pplane5(''eqpt'')','separator','on');
         
         menu(4) = uimenu(solmenu,...
            'label','List computed equilibrium points.',...
            'call','pplane5(''eqptlist'')');
         
         mestunst= uimenu(solmenu,'label','Plot stable and unstable orbits.',...
            'call','pplane5(''stunst'')','interrupt','on');
         
         nullcall = ['me = gcbo;',...
               'dud = get(gcf,''user'');',...
               'handx = [dud.arr.hx; dud.arr.hy];',...	
               'handr = [dud.arr.hr; dud.arr.hth];',...	
               'lab = get(me,''label'');',...
               'switch lab,',...
               '  case ''Plot nullclines.'',',...
               '    set(handx,''vis'',''on'');',...
               '    set(handr,''vis'',''off'');',...
               '    set(me,''label'',''Delete nullclines.'');',...
               '    set(findobj(''tag'',''rnull''),''label'',''Plot polar nullclines.'');',...
               '  case ''Delete nullclines.'','...
               '    set(handx,''vis'',''off'');',...
               '    set(me,''label'',''Plot nullclines.'');',...
               'end'];
         
         rnullcall = ['me = gcbo;',...
               'dud = get(gcf,''user'');',...
               'handx = [dud.arr.hx; dud.arr.hy];',...	    
               'handr = [dud.arr.hr; dud.arr.hth];',...	
               'lab = get(me,''label'');',...
               'switch lab,',...
               '  case ''Plot polar nullclines.'',',...
               '    set(handr,''vis'',''on'');',...
               '    set(handx,''vis'',''off'');',...
               '    set(me,''label'',''Delete polar nullclines.'');',...
               '    set(findobj(''tag'',''null''),''label'',''Plot nullclines.'');',...
               '  case ''Delete polar nullclines.'','...
               '    set(handr,''vis'',''off'');',...
               '    set(me,''label'',''Plot polar nullclines.'');',...
               'end'];
         
         menunull = uimenu(solmenu,'label','Plot nullclines.',...
            'call',nullcall,'separator','on','tag','null');
         
         if dud.egg
            menurnull = uimenu(solmenu,'label','Plot polar nullclines.',...
               'call',rnullcall,'separator','off','tag','rnull');
            metest = uimenu(solmenu,'label','Test case',...
               'call','pplane5(''test case'')',...
               'separator','on','tag','testcase');
            
         end      
         
         medall = uimenu(emenu,'label','Erase all solutions.',...
            'call','pplane5(''dall'')','separator','on');
         
         medel = uimenu(emenu,'label','Delete a graphics object.',...
            'call','pplane5(''delete'')','visible','on');
         
         menutext = uimenu(emenu,'label','Enter text on the Display Window.',...
            'call','pplane5(''text'')');
         
         meset   = uimenu(menu(1),'label','Settings.',...
            'call','pplane5(''settings'')');
         
         mesolve = uimenu(menu(1),'label','Solver.');
         
         solverstr = ['ud = get(gcf,''user'');',...
               'me = gcbo;',...
               'meud = get(me,''user'');',...
               'ud.settings.refine = meud.refine;',...
               'ud.settings.tol = meud.tol;',...
               'ud.settings.solver = meud.solver;',...
	       'ud.settings.stepsize = meud.stepsize;',...
               'set(ud.solver,''checked'',''off'');',...
               'set(me,''checked'',''on'');',...
               'set(gcf,''user'',ud);',...
               'ppset = findobj(''name'',''PPLANE5 Setup'');',...
               'sud = get(ppset,''user'');',...
               'sud.settings = ud.settings;',...
               'set(ppset,''user'',sud);',...
               'pplane5(''settings'');'];
         
      solver = dud.settings.solver;
      dpset.refine = 8;
      dpset.tol = dud.settings.tol;
      dpset.solver = 'Dormand Prince';
      dpset.stepsize = dud.settings.stepsize;
      dpset.hmax = 0;
      if strcmp(solver,'Dormand Prince')
	dpch = 'on';
      else
	dpch = 'off';
      end
      rk4set.refine = 1;
      rk4set.tol = dud.settings.tol;
      rk4set.solver = 'Runge-Kutta 4';
      rk4set.stepsize = dud.settings.stepsize;
      rk4set.hmax = 0;
      if strcmp(solver,'Runge-Kutta 4')
	rk4ch = 'on';
      else
	rk4ch = 'off';
      end
      ode15sset.refine = 1;
      ode15sset.tol = dud.settings.tol;
      ode15sset.solver = 'ode15s';
      ode15sset.stepsize = dud.settings.stepsize;
      ode15sset.hmax = 0;
      if strcmp(solver,'ode15s')
	ode15sch = 'on';
      else
	ode15sch = 'off';
      end
      ode23sset.refine = 1;
      ode23sset.tol = dud.settings.tol;
      ode23sset.solver = 'ode23s';
      ode23sset.stepsize = dud.settings.stepsize;
      ode23sset.hmax = 0;
      if strcmp(solver,'ode23s')
	ode23sch = 'on';
      else
	ode23sch = 'off';
      end
      ode113set.refine = 1;
      ode113set.tol = dud.settings.tol;
      ode113set.solver = 'ode113';
      ode113set.stepsize = dud.settings.stepsize;
      ode113set.hmax = 0;
      if strcmp(solver,'ode113')
	ode113ch = 'on';
      else
	ode113ch = 'off';
      end
      ode23set.refine = 1;
      ode23set.tol = dud.settings.tol;
      ode23set.solver = 'ode23';
      ode23set.stepsize = dud.settings.stepsize;
      ode23set.hmax = 0;
      if strcmp(solver,'ode23')
	ode23ch = 'on';
      else
	ode23ch = 'off';
      end
      ode45set.refine = 8;
      ode45set.tol = dud.settings.tol;
      ode45set.solver = 'ode45';
      ode45set.stepsize = dud.settings.stepsize;
      ode45set.hmax = 0;
      if strcmp(solver,'ode45')
	ode45ch = 'on';
      else
	ode45ch = 'off';
      end
         dud.solver(1) = uimenu(mesolve,'label','Dormand Prince',...
            'checked',dpch,...
            'call',solverstr,'user',dpset);
         
         dud.solver(2) = uimenu(mesolve,'label','ode45',...
            'checked',ode45ch,...
            'call',solverstr,'user',ode45set);
         
         dud.solver(3) = uimenu(mesolve,'label','ode23',...
            'checked',ode23ch,...
            'call',solverstr,'user',ode23set);
         
         dud.solver(4) = uimenu(mesolve,'label','ode113',...
            'checked',ode113ch,...
            'call',solverstr,'user',ode113set);
         
         dud.solver(5) = uimenu(mesolve,'label','ode15s',...
            'checked',ode15sch,...
            'call',solverstr,'user',ode15sset);
         
         dud.solver(6) = uimenu(mesolve,'label','ode23s',...
            'checked',ode23sch,...
            'call',solverstr,'user',ode23sset);
         
         dud.solver(7) = uimenu(mesolve,'label','Runge-Kutta 4',...
            'checked',rk4ch,...
            'call',solverstr,'user',rk4set);
	
         if ~dud.egg
            set(dud.solver(2:5),'vis','off');
         end
         
         plotch = [
            'dud = get(gcf,''user'');',...
               'aud = get(dud.axes,''user'');',...
               'if aud.plot,',...
               '   aud.plot = 0;',...
               '   set(gcbo,''label'',''Plot while computing'');',...
               'else,',...
               '   aud.plot = 1;',...
               '   set(gcbo,''label'',''Do not plot while computing'');',...
               'end,',...
               'set(dud.axes,''user'',aud);'];
         
         medir = uimenu(menu(1),'label','Solution direction.');
         
         directionstr = ['ud = get(gcf,''user'');',...
               'me = gcbo;',...
               'ud.dir = get(me,''user'');',...
               'set(ud.direction,''checked'',''off'');',...
               'set(me,''checked'',''on'');',...
               'set(gcf,''user'',ud);'];
         
         dud.direction(1) = uimenu(medir,'label','Both',...
            'checked','on',...
            'user',0,...
            'call',directionstr);
         dud.dir = 0;
         
         dud.direction(2) = uimenu(medir,'label','Forward',...
            'user',1,...
            'call',directionstr);
         
         dud.direction(3) = uimenu(medir,'label','Back',...
            'user',-1,...
            'call',directionstr);
         
         meplot = uimenu(menu(1),'label','Do not plot while computing',...
            'call',plotch,'separator','on');
         
         menu(6) = uimenu(menu(1),'label','Make the Display Window inactive.',...
            'call','pplane5(''hotcold'')','separator','on');
         
         dud.menu = menu;
         
         if exist('RICE') == 2
            ricemenu = uimenu (menu(1),'label','Set default printer.',...
               'call','pplane5 (''riceprint'')','separator','on');
         end
         
         
         set(gcf,'WindowButtonDownFcn','pplane5(''down'')');
         hh1 = [dud.axes,dud.title.axes,nframe,dud.notice,dbutt([1 2 3])];
         set(hh1,'units','norm');
         hh2 = [nframe,dud.notice,dbutt(2:3),dud.axes,dud.menu(1),meplot,megraph];
         set(hh2,'visible','on');
         
         set(ppdisp,'vis','on');
         dud.printstr = 'print -noui';
      end   % if (~isempty(ppdisp))  & else
      ppdispa = dud.axes;
      axes(ppdispa);
      cla
      xlabel(txname);
      ylabel(tyname);
      % Initialize the window matrix.
      
      set(findobj('tag','zbmenu'),'enable','off');
    tstr1 = [txname,' '' = ', txderstr];
    tstr2 = [tyname,' '' = ', tyderstr];
    tstr = str2mat(tstr1,tstr2);
    dud.tstr = tstr;
      
      
      
      k = find(pflag);
      if ~isempty(k)
         lk = length(k);
         switch lk
         case 1
            pstr1 = [pstring(k);{' '}];
         case 2
            pstr1 = pstring(k);
         case 3
            pstr1 = pstring(k([2,3]));
            pstr2 = [pstring(k(1));{' '}];
         case 4
            pstr1 = pstring(k([2,4]));
            pstr2 = pstring(k([1,3]));
         end
      end
      set(dud.title.eq,'string',tstr);
      set(dud.title.p1,'string',pstr1);
      ext = get(dud.title.p1,'extent');
      pos = get(dud.title.p1,'pos');
      p1 = min(.9, .93 - ext(3));
      pos(1) = p1;
      set(dud.title.p1,'pos',pos);
      set(dud.title.p2,'string',pstr2);	
      ext = get(dud.title.p2,'extent');
      pos = get(dud.title.p2,'pos');
      pos(1) = min(.8, p1 - ext(3)-0.02);
      set(dud.title.p2,'pos',pos);
      
      % Initialize important information as user data.
      
      dud.function = dfcn;
      dud.solhand = [];	% Handles to solution curves and equilibrium points.
      dud.arr = [];	% Handles for the direction and vector fields.
      dud.eqpts = [];	% Equilibrium point data.
      dud.wmat = [];
      dud.color = sud.color;
      set(ppdisp,'user',dud);
      if strcmp(dud.syst.fieldtype,'nullclines')
         menunull = findobj('tag','null');      
         set(menunull,'enable','off');
      end      
      ud.y = zeros(2,1);
      ud.i = 0;
      ud.line = 0;
      wind = dud.syst.wind(:);
      dwind = [wind(1); wind(3); -wind(2); -wind(4)];
      ud.DY = [wind(2)-wind(1); wind(4)-wind(3)];
      ud.cwind = dwind - dud.settings.magn*[ud.DY;ud.DY];
      ud.R = zeros(2,2);
      ud.rr = zeros(2,2);
      ud.perpeps = 0;
      ud.paraeps = 0;
      ud.sinkeps = 0;
      ud.turn = zeros(2,10);
      ud.tk = 0;
      ud.stop = 0;
      ud.gstop = 1;
      ud.plot = 1;
      set(dud.axes,'user',ud);
      tc = findobj('tag','testcase');
      if ~isempty(tc)    
         if strcmp(dud.syst.name,'default system')
            set(tc,'vis','on');
         else
            set(tc,'vis','off');
         end
      end    
      ppkbd = findobj('name','PPLANE5 Keyboard input','vis','on');
      if ~isempty(ppkbd),pplane5('kbd'),end
      set(dud.title.axes,'handle','off');  
      pplane5('dirfield',ppdisp);		
   end
   
elseif strcmp(action,'linear')
   
   % Initialize linearization window.
   
   % Linear takes the information from the EQPT window and initializes
   % the Linear Display window if it already exists.  If the Linear 
   % Display window does not exist, it builds one.
   
   ppeqpt = gcf;
   
   % Get the information from PPLANE5 Equilibrium ...
   
   sud = get(ppeqpt,'user');
   WINvect = [-1 1 -1 1];
   jac = sud.jac;
   type = sud.type;
   vectors = sud.vectors;
   ppdisp = findobj('name','PPLANE5 Display');
   pdud = get(ppdisp,'user');
   settings = pdud.settings;
   system = sud.system;
   
   % Find PPLANE5 Linearization if it exists.	
   % If PPLANE5 Linearization exists, update it.  If it does not build it.
   
   pplin= findobj('name','PPLANE5 Linearization');
   if (~isempty(pplin))
      figure(pplin);
      dud = get(pplin,'user');
      dud.syst = system;
      dud.settings = settings;
      dfcn = dud.function;
   else
      pplin = figure('name','PPLANE5 Linearization',...
         'numb','off',...
         'interrupt','on',...
         'visible','off',...
         'tag','PPLANE5');
      pplane5('figdefault',pplin);
      dud = get(pplin,'user');
      dud.dir = 0;
      dud.syst = system;
      dud.settings = settings;
      dfcn = '';
      fs = dud.fontsize;
      ssize = dud.ssize;
      r = ssize/10;
      ppaxw = 300;	% Default axes width
      ppaxh = 300;	% Default axes height
      ppaxl = 45;	% Default axes left
      buttw = 40;     % Default button width
      titleh = 45;    % Default title height
      ppaxb = 35;
      ppdh = ppaxb + ppaxh + titleh;
      uni = get(0,'units');
      set(0,'units','pixels');
      ss = get(0,'screensize');
      set(0,'units',uni);
      sw = ss(3);sh = ss(4);
      if r*ppdh > sh -80;
         r = (sh-80)/ppdh;
         fs = fs*r;
         lw = 0.5*r;
         set(gcf,'defaultaxesfontsize',fs,'defaultaxeslinewidth',lw);
         set(gcf,'defaulttextfontsize',fs);
         set(gcf,'defaultlinelinewidth',lw);
         set(gcf,'defaultuicontrolfontsize',fs*0.9);
      end
      
      axpos = r*[ppaxl,ppaxb,ppaxw,ppaxh];
      
      % Set up the plot axes.
      
      dud.axes = axes('units','pix',...
         'position',axpos,...
         'next','add',...
         'box','on',...
         'interrupt','on',...
         'xgrid','on',...
         'ygrid','on',...
         'drawmode','fast',...
         'visible','off',...
         'tag','linear axes');
      
      % Set up the title.
      
      dud.title.axes = axes('box','off','xlim',[0 1],'ylim',[0 1],...
         'units','pix','vis','off','xtick',[-1],'ytick',[-1],...
         'xticklabel','','yticklabel','');
      
      dud.title.eq = text(0.07,0.5,' ','vert','middle');
      dud.title.p1 = text(0.75,0.5,' ','vert','middle');
      dud.title.p2 = text(0.65,0.5,' ','vert','middle');
      tstr = {'x_2';'y^2'};
      set(dud.title.eq,'string',tstr,'units','pix');
      ext = get(dud.title.eq ,'extent');
      titleh = ext(4)+15*r;
      set(dud.title.eq,'units','data');
      taxpos = r*[ppaxl,ppaxb+ppaxh,ppaxw,titleh];
      set(dud.title.axes,'pos',taxpos,'color',get(gcf,'color'));
      
      % Finish the positions.
      
      % 		buttw = r*buttw;
      % 		butth = fs+10*r;
      % 		buttl = ppaxl + ppaxw + 5;
      ppdw = r*(ppaxl + ppaxw + 35);
      ppdh = r*(ppaxb+ppaxh)+titleh;
      
      ppdleft = sw-ppdw-60;
      ppdbot = 60;
      ppdpos = [ppdleft,ppdbot,ppdw,ppdh];
      set(pplin,'resize','on');
      set(pplin,'pos',ppdpos);
      
      Arrflag = system.fieldtype;
      
      axpos = r*[ppaxl,ppaxb,ppaxw,ppaxh];
      
      stopstr = 'aud = get(gca,''user'');aud.stop = 4;set(gca,''user'',aud);';
      
      stoppos = [axpos(1)+axpos(3)-r*10, axpos(2)+axpos(4)+r*10, r*40,fs+10*r];    
      
      dud.butt  = uicontrol('style','push',...
         'pos',stoppos,...
         'string','Stop','call',stopstr,'vis','off','tag','stop');
      
      fmenu = findobj(gcf,'label','File');  % uimenu('label','File');
      emenu = findobj(gcf,'label','Edit');  % uimenu('label','Edit');
      
      solmenu = uimenu('label','Solutions');
      
      optmenu = uimenu('label','Options','visible','off');
      
      meprint = uimenu(fmenu,'label','Print PPLANE5 Linearization',...
         'HandleVis','on', ...
         'Accelerator','P', ...
         'Callback','print');
      
      meprsetup =   uimenu(fmenu,'label','Page Position ...',...
         'Handlevis','on', ...
         'Callback','pagedlg(gcf)'); 
      
      meprdlg =   uimenu(fmenu,'label','Print ...',...
         'Handlevis','on', ...
         'Callback','printdlg(gcf)');  
      
      meclose = uimenu(fmenu,'label','Close',...
         'Handlevisibility','on', ...
         'Accelerator','W', ...
         'Callback','close');
      
      mequit = uimenu(fmenu,'label','Quit PPLANE5',...
         'call','pplane5(''quit'')','separator','on');
      
      menukey = uimenu(solmenu,'label','Keyboard input.','call',...
         ['pplane5(''kbd'',',num2str(dud.axes),')'],'vis','off');
      
      mesev   = uimenu(solmenu,'label','Plot several solutions.',...
         'call','pplane5(''several'')');
      
      
      fundcall = ['dud = get(gcf,''user'');',...      
            'col = dud.color.sep;',...
            'vect = dud.vectors;',...
            'h(1) = plot(2*[vect(1,1),-vect(1,1)],2*[vect(2,1),-vect(2,1)]);',...
            'h(2) = plot(2*[vect(1,2),-vect(1,2)],2*[vect(2,2),-vect(2,2)]);',...
            'set(h,''color'',col);',...	    
            'dud.solhand = [dud.solhand;h(:)];'...
            'set(gcf,''user'',dud);'];
      
      dud.menu = uimenu(solmenu,'call',fundcall);
      
      medall = uimenu(emenu,'label','Erase all solutions.',...
         'call','pplane5(''dall'')');
      
      medel = uimenu(emenu,'label','Delete a graphics object.',...
         'call','pplane5(''delete'')','visible','on');
      
      menutext = uimenu(emenu,'label','Enter text on the Display Window.',...
         'call','pplane5(''text'')');
      
      mehc = uimenu(optmenu,'label','Make the Display Window inactive.',...
         'call','pplane5(''hotcold'')','separator','on');
      
      set(pplin,'WindowButtonDownFcn',['pplane5(''down'',',num2str(pplin),')']);
      hh1 = [dud.axes,dud.title.axes];
      set(hh1,'units','norm');
      set(dud.axes,'visible','on');
      set(pplin,'vis','on');
      dud.printstr = 'print';
   end   % if (~isempty(pplin))  & else
   
   switch type
   case 1
      vstr = 'on';
      lstr = 'Plot stable and unstable orbits.';
   case {2,3}
      vstr = 'on';
      lstr = 'Plot fundamental modes.';
   otherwise
      vstr = 'off';
      lstr='';
   end
   set(dud.menu,'label',lstr,'vis',vstr);
   
   pplina = dud.axes;
   dud.vectors = sud.vectors;
   axes(pplina);
   cla
   xlabel('u');
   ylabel('v');
   
   % The title strings.
   
   tstr = {'u'' = A u + B v';'v'' = C u + D v'};
   dud.tstr = tstr;
   pstr2 = {['A = ',num2str(jac(1,1))];['C = ',num2str(jac(2,1))]};
   pstr1 = {['B = ',num2str(jac(1,2))];['D = ',num2str(jac(2,2))]};
   
   set(dud.title.eq,'string',tstr);
   set(dud.title.p1,'string',pstr1);
   ext = get(dud.title.p1,'extent');
   pos = get(dud.title.p1,'pos');
   p1 = min(.9, .93 - ext(3));
   pos(1) = p1;
   set(dud.title.p1,'pos',pos);
   set(dud.title.p2,'string',pstr2);	
   ext = get(dud.title.p2,'extent');
   pos = get(dud.title.p2,'pos');
   pos(1) = min(.8, p1 - ext(3)-0.02);
   set(dud.title.p2,'pos',pos);
   
   if (~strcmp(dfcn,'') & exist(dfcn)==2) delete([dfcn,'.m']);end
   tee = clock;
   tee = ceil(tee(6)*100);
   dfcn=['pptp',num2str(tee)];
   fcnstr = ['function ypr = ',dfcn,'(t,y)\n\n'];
   commstr = '%%%% Created by PPLANE5\n\n';
   astr = ['A = [',num2str(jac(1,1)),',',num2str(jac(1,2)),';',...
         num2str(jac(2,1)),',',num2str(jac(2,2)),'];\n'];
   dstr = 'ypr = A * y;';
   ppf = fopen([dfcn,'.m'],'w');
   fprintf(ppf,fcnstr);
   fprintf(ppf,commstr);
   fprintf(ppf,astr);
   fprintf(ppf,dstr);
   fclose(ppf);
   
   % Initialize important information as user data.
   
   dud.function = dfcn;
   dud.solhand = [];	% Handles to solution curves and equilibrium points.
   dud.arr = [];		% Handles for the direction and vector
   % fields.
   dud.eqpts = [];		% Equilibrium point data.
   dud.notice = 0;
   dud.syst.wind = [-1 1 -1 1]';
   dud.color = sud.color;
   set(pplin,'user',dud);
   
   ud.y = zeros(2,1);
   ud.i = 0;
   ud.line = 0;
   wind = [-1 1 -1 1]';
   dwind = [wind(1); wind(3); -wind(2); -wind(4)];
   ud.DY = [wind(2)-wind(1); wind(4)-wind(3)];
   ud.cwind = dwind - dud.settings.magn*[ud.DY;ud.DY];
   ud.R = zeros(2,2);
   ud.rr = zeros(2,2);
   ud.perpeps = 0;
   ud.paraeps = 0;
   ud.sinkeps = 0;
   ud.turn = zeros(2,10);
   ud.tk = 0;
   ud.stop = 0;
   ud.gstop = 1;
   ud.plot = 1;
   set(dud.axes,'user',ud);
   ppkbd = findobj('name','PPLANE5 Keyboard input','vis','on');
   if ~isempty(ppkbd),pplane5('kbd'),end
   pplane5('dirfield',pplin);		
   plot(0,0,...
      'markerfacecolor',dud.color.eqpt,...
      'markeredgecolor',dud.color.eqpt,...
      'marker','o');
   
elseif strcmp(action,'dirfield')
   
   % 'dirfield' computes and plots the field elements.  This is the entry
   % point both from 'display' and from later commands that require the
   % recomputation of the field elements.
   
   % Find PPLANE5 Display and get the user data.
   
   disph = input1;   % This could be ppdisp or pplin.
   
   dud = get(disph,'user');
   color = dud.color;
   
   dfcn = dud.function;
   ppdispa = dud.axes;
   WINvect = dud.syst.wind;
   settings = dud.settings;
   notice = dud.notice;
   if notice
      nstr = get(notice,'string');
      nstr(1:4)=nstr(2:5);
      nstr{5,1} = 'Computing the field elements.';
      set(notice,'string',nstr);
      
      % Augment the window matrix
      
      wmat = dud.wmat;
      wrows = size(wmat,1);
      wflag = 0;
      for k = 1:wrows
         if wmat(k,:)==WINvect
            wflag = 1;
         end
      end
      if wflag == 0
         wmat = [wmat;WINvect];
         dud.wmat = wmat;
      end
      if wrows 
         hhh = findobj('tag','zbmenu');
         set(hhh,'enable','on');
      end
   end
   
   
   Xmin = WINvect(1);
   Xmax = WINvect(2);
   Ymin = WINvect(3);
   Ymax = WINvect(4);
   
   N = dud.syst.npts;
   if strcmp(get(disph,'name'),'PPLANE5 Linearization') 
      N = floor(3*N/4);
   end    
   deltax=(Xmax - Xmin)/(N-1);
   deltay=(Ymax - Ymin)/(N-1);
   
   % Set up the display window.
   
   Dxint=[Xmin-deltax,Xmax+deltax];
   Dyint=[Ymin-deltay,Ymax+deltay];
   
   % Set up the original mesh.
   
   XXXg=Xmin + deltax*[0:N-1];
   YYYg=Ymin + deltay*[0:N-1];
   
   [Xx,Yy]=meshgrid(XXXg,YYYg);
   
   % Calculate the line and vector fields.
   
   Xx=Xx(:);Yy=Yy(:);
   Ww = zeros(size(Xx));
   Ww = feval(dfcn,0,[Xx';Yy']);
   Vv = Ww(1,:) + Ww(2,:)*sqrt(-1);
   Vv = Vv.';
   Arrflag = dud.syst.fieldtype;
   
   
   mgrid = Xx+Yy.*sqrt(-1); % mgrid = mgrid(:);
   zz=Vv.';
   sc = min(deltax,deltay);
   
   arrow=[-1,1].';
   zzz=sign(zz);
   scale = sqrt((real(zzz)/deltax).^2+(imag(zzz)/deltay).^2);
   ww = (zzz == 0);
   scale = scale + ww;
   aa1 = 0.3*arrow*(zzz./scale)+ones(size(arrow))*(mgrid.');
   [r,c] = size(aa1);
   aa1=[aa1;NaN*ones(1,c)];
   aa1=aa1(:);
   
   arrow = [0,1,.7,1,.7].' + [0,0,.25,0,-.25].' * sqrt(-1);
   zz=sign(zz).*((abs(zz)).^(1/3));
   scale = 0.9*sc./max(max(abs(zz)));
   aa2 = scale*arrow*zz +ones(size(arrow))*(mgrid.');
   [r,c] = size(aa2);
   aa2=[aa2;NaN*ones(1,c)];
   aa2=aa2(:);
   axes(ppdispa);
   
   arr = dud.arr;	% Delete the old field data.
   if isstruct(arr)
      hand = [arr.hx;arr.hy;arr.lines;arr.arrows;arr.barrows];  
      delete(hand);
   end
   NN = N;
   N = 50; k = ceil(N/NN);
   deltax=(Xmax - Xmin)/(N-1);
   deltay=(Ymax - Ymin)/(N-1);
   
   % Set up the original mesh.
   
   XXXg=Xmin + deltax*[-k:N+k];
   YYYg=Ymin + deltay*[-k:N+k];
   
   [Xx,Yy]=meshgrid(XXXg,YYYg);
   
   % Calculate the line and vector fields.
   
   Xxx=Xx(:);Yyy=Yy(:);
   Ww = zeros(size(Xxx));
   Ww = feval(dfcn,0,[Xxx';Yyy']);
   DX = Ww(1,:)';
   DY = Ww(2,:)';
   minx = min(DX);
   maxx = max(DX);
   miny = min(DY);
   maxy = max(DY);
   DR = Xxx.*DX + Yyy.*DY;
   DTheta = Xxx.*DY - Yyy.*DX;
   minr = min(DR);
   maxr = max(DR);
   minth = min(DTheta);
   maxth = max(DTheta);
   DX = reshape(DX,N+2*k+1,N+2*k+1);
   DY = reshape(DY,N+2*k+1,N+2*k+1);
   DR = reshape(DR,N+2*k+1,N+2*k+1);
   DTheta = reshape(DTheta,N+2*k+1,N+2*k+1);
   if minx < 0 & 0 < maxx
      [Cx,hx] = contour(Xx,Yy,DX,[0,0],'--');
      set(hx,'visible','off','color',color.xcline,'linestyle','--');
   else
      hx = zeros(0,1);
   end
   if miny < 0 & 0 < maxy
      [Cy,hy] = contour(Xx,Yy,DY,[0,0],'--');
      set(hy,'visible','off','color',color.ycline,'linestyle','--');
   else
      hy = zeros(0,1);
   end
   if minr < 0 & 0 < maxr
      [Cr,hr] = contour(Xx,Yy,DR,[0,0],'--');
      set(hr,'visible','off','color',color.xcline,'linestyle','--');
   else
      hr = zeros(0,1);
   end
   if minth < 0 & 0 < maxth
      [Cth,hth] = contour(Xx,Yy,DTheta,[0,0],'--');
      set(hth,'visible','off','color',color.ycline,'linestyle','--');
   else
      hth = zeros(0,1);
   end
   arrh1 = plot(real(aa1),imag(aa1),'color',color.arrows,'visible','off');
   arrh2 = plot(real(aa2),imag(aa2),'color',color.arrows,'visible','off');
   
   % We plot both the line field and the vector field.  Then we
   % control which is seen by manipulating the visibility.
   Xmin = Xmin - k*deltax;
   Xmax = Xmax + k*deltax; 
   Ymin = Ymin - k*deltay;
   Ymax = Ymax + k*deltay;  
   Zz = (Xx-Xmin).*(Xmax-Xx).*(Yy-Ymin).*(Ymax-Yy).*abs(DX).*abs(DY);
   Zzc = Zz(2:N+k+2,2:N+k+2);
   Xxc = Xx(2:N+k+2,2:N+k+2);
   Yyc = Yy(2:N+k+2,2:N+k+2);
   DXc = DX(2:N+k+2,2:N+k+2);
   DYc = DY(2:N+k+2,2:N+k+2);
   Zzl = Zz(1:N+k+1,2:N+k+2);
   Zzr = Zz(3:N+k+3,2:N+k+2);
   Zzd = Zz(2:N+k+2,1:N+k+1);
   Zzu = Zz(2:N+k+2,3:N+k+3);
   Zzdl = Zz(1:N+k+1,1:N+k+1);
   Zzdr = Zz(3:N+k+3,1:N+k+1);
   Zzul = Zz(1:N+k+1,3:N+k+3);
   Zzur = Zz(3:N+k+3,3:N+k+3);
   kk = find((Zzc>Zzl)&(Zzc>Zzr)&(Zzc>Zzd)&(Zzc>Zzu)...
      &(Zzc>Zzdl)&(Zzc>Zzdr)&(Zzc>Zzul)&(Zzc>Zzur));
   Xxx = Xxc(kk);
   Yyy = Yyc(kk);
   DXx = DXc(kk);
   DYy = DYc(kk);
   ll = length(hx);
   Xxx = Xxx(:);
   Yyy = Yyy(:);
   DXx = DXx(:);
   DYy = DYy(:);
   if ~exist('both')
      Xxx = zeros(0,1);
      Yyy = zeros(0,1);
      DXx = zeros(0,1);
      DYy = zeros(0,1);
   end
   for j = 1:ll
      xd = get(hx(j),'xdata');xd = xd(:);
      yd = get(hx(j),'ydata');yd = yd(:);
      zxd = feval(dfcn,0,[xd';yd']);
      yxd = zxd(2,:);yxd = yxd(:);
      ayxd = abs(yxd).*(xd-Xmin).*(Xmax-xd).*(yd-Ymin).*(Ymax-yd);
      NNN = length(ayxd);
      ayxdc = ayxd(2:NNN-1);
      yxdc = yxd(2:NNN-1);
      ayxdl = ayxd(1:NNN-2);
      ayxdr = ayxd(3:NNN);
      lll = find((ayxdc>ayxdl)&(ayxdc>ayxdr));
      Xx = xd(lll+1);
      Yy = yd(lll+1);
      Xxx = [Xxx;Xx(:)];
      Yyy = [Yyy;Yy(:)];
      yxdcp = yxdc(lll);
      yxdcp = yxdcp(:);
      DXx = [DXx;zeros(size(yxdcp))];
      DYy = [DYy;yxdcp];
   end
   ll = length(hy);
   for j = 1:ll
      xd = get(hy(j),'xdata');xd = xd(:);
      yd = get(hy(j),'ydata');yd = yd(:);
      zxd = feval(dfcn,0,[xd';yd']);
      yxd = zxd(1,:);yxd = yxd(:);
      ayxd = abs(yxd).*(xd-Xmin).*(Xmax-xd).*(yd-Ymin).*(Ymax-yd);
      NNN = length(ayxd);
      ayxdc = ayxd(2:NNN-1);
      yxdc = yxd(2:NNN-1);
      ayxdl = ayxd(1:NNN-2);
      ayxdr = ayxd(3:NNN);
      lll = find((ayxdc>ayxdl)&(ayxdc>ayxdr));
      Xx = xd(lll+1);
      Yy = yd(lll+1);
      Xxx = [Xxx;Xx(:)];
      Yyy = [Yyy;Yy(:)];
      yxdcp = yxdc(lll);
      yxdcp = yxdcp(:);
      DXx = [DXx;yxdcp];
      DYy = [DYy;zeros(size(yxdcp))];
   end
   mgrid = Xxx(:) + Yyy(:)*sqrt(-1);
   Vv = DXx(:) + DYy(:)*sqrt(-1);
   Vv = sign(Vv);       %.*((abs(Vv)).^(1/3));
   sc = (N/20)*min(deltax,deltay);
   aa3 = sc*arrow*(Vv.') +ones(size(arrow))*(mgrid.');
   [r,c] = size(aa3);
   aa3 = [aa3;NaN*ones(1,c)];
   aa3=aa3(:);
   arrh3 = plot(real(aa3),imag(aa3),'color',color.narrows,'visible','off');
   
   switch Arrflag
   case 'nullclines'
      set([hx;hy;arrh3],'vis','on');
   case 'lines'
      set(arrh1,'visible','on');
   case 'arrows'
      set(arrh2,'visible','on');
   end
   dud.arr.lines = arrh1;	   % Save the handles for later use.
   dud.arr.arrows = arrh2;	   % Save the handles for later use.
   dud.arr.hx = hx;
   dud.arr.hy = hy;
   dud.arr.barrows = arrh3;
   dud.arr.hr = hr;
   dud.arr.hth = hth;
   
   menull = findobj('tag','null');
   mernull = findobj('tag','rnull');
   if strcmp(get(menull,'label'),'Delete nullclines.')
      set([hx;hy],'vis','on');
   elseif strcmp(get(mernull,'label'),'Delete polar nullclines.')
      set([hr;hth],'vis','on');
   end
   if notice
      nstr = get(notice,'string');
      nstr(1:4) = nstr(2:5);
      nstr{5,1} = 'Ready.';
      set(notice,'string',nstr);
   end
   set(disph,'user',dud);
   axis([Dxint,Dyint]);
   
elseif strcmp (action,'riceprint')
   
   printcode = [	
      'print             ';
      'print -Pbaker     ';
      'print -Pbrown     ';
      'print -Phanszen   ';
      'print -Pjones     ';
      'print -Plovett    ';
      'print -Psid       ';
      'print -Pwiess     ';
      'print -Pwrc       ';
      'print -Pfondren1ps';
      'print -Pfondren2ps';
      'print -Pmd105     '
   ];
   
   printname = [
      'The printer in this location';
      'Baker College               ';
      'Brown College               ';
      'Hanszen College             ';
      'Jones College               ';
      'Lovett College              ';
      'Sid Rich College            ';
      'Wiess College               ';
      'Will Rice College           ';
      'Fondren (front)             ';
      'Fondren (back)              ';
      'Mudd 105                    ';
   ];
   
   NN = size(printname,1);
   
   dud = get(gcf,'user');
   
   txtl = 150;txth = 21;
   
   pbwidth = 80;
   framewidth = txtl+pbwidth+20;
   frht = 10 + NN*txth +2*txth +10;
   frleft = 10;
   frbot = 20 + txth;
   frpos = [frleft frbot framewidth frht];
   figwidth = framewidth + 20;
   
   dfrice = figure('pos',[150, 50, figwidth, frbot+frht + 10],...
      'name','PPLANE5 Rice printer selection','numb','off',...
      'vis','off','tag','PPLANE5');
   
   pplane5('figdefault',dfrice);
   
   set(dfrice,'user',dud.printstr);		
   
   
   riceframe = uicontrol('style','frame','pos',frpos,'vis','off');
   printstr = dud.printstr;
   maxval = 1:NN;
   
   for j=1:NN
      
      rad(j) = uicontrol('style','radio',...
         'pos',[(frleft+8), frbot+10+(NN-j)*txth pbwidth+txtl txth],...
         'string',printname(j,:),'max',j,...
         'value',strcmp(printstr,deblank(printcode(j,:)))*maxval(j),...
         'tag',deblank(printcode(j,:)),'hor','left','vis','off');
   end    
   
   
   for i=1:NN
      set(rad(i),'user',rad(:,[1:(i-1),(i+1):NN]));
   end
   
   callrad = [
      'ppdisp = findobj(''name'',''PPLANE5 Display'');'...
         'dud = get(ppdisp,''user'');',...
         'me = get(gcf,''currentobject'');',...
         'kk = get(me,''max'');',...
         'set(get(me,''user''),''value'',0),',...
         'set(me,''value'',kk),',...
         'dud.printstr = get(me,''tag'');',...
         'set(ppdisp,''user'',dud);'];
   
   set(rad,'call',callrad);
   textpos = [frleft+8, frbot+10+(NN+0.5)*txth,pbwidth+txtl,txth];
   textstr = 'Select the printer of your choice';
   textw=uicontrol('style','text','pos',textpos,'string',textstr,...
      'hor','left','vis','off');
   bwidth=80;
   bpos = [(figwidth-bwidth)/2,10,bwidth,txth];
   offbutt = uicontrol('style','push','pos',bpos,'string','Go away',...
      'call','delete(gcf)');
   
   set(findobj(gcf,'type','uicontrol'),'units','normal');
   set(gcf,'resize','on');	
   set(gcf,'vis','on');
   set(get(gcf,'child'),'vis','on');
   
elseif strcmp(action,'hotcold')
   
   % 'hotcold' is the callback for the menu selection that makes the
   % Display Window active or inactive.
   
   ppdisp = gcf;
   dud = get(ppdisp,'user');
   nstr = get(dud.notice,'string');
   nstr(1:4) = nstr(2:5);
   mehc = dud.menu(6);
   if (findstr(get(mehc,'label'),'inactive'))
      set(ppdisp,'WindowButtonDownFcn',' ');
      set(mehc,'label','Make the Display Window active.');
      nstr{5,1} = 'The Display Window is not active.';
      set(dud.notice,'string',nstr);
   else
      set(ppdisp,'WindowButtonDownFcn','pplane5(''down'')');
      set(mehc,'label','Make the Display Window inactive.');
      nstr{5,1} = 'The Display Window is active.';
      set(dud.notice,'string',nstr);
   end
   
elseif strcmp(action,'down')
   
   % 'down' is the Window Button Down call.  It starts the computation of
   % solutions from a click of the mouse.
   
   disph = gcf;
   seltype = get(disph,'selectiontype');
   if strcmp(seltype,'alt')
      pplane5('zoom');
      return
   elseif strcmp(seltype,'extend')
      pplane5('zoomsqd');
      return
   end
   dud = get(disph,'user');
   ax = dud.axes;
   ch = findobj('type','uicontrol','enable','on');
   %  mh = findobj('type','uimenu','enable','on');
   set(ch,'enable','inactive');
   %  set(mh,'enable','off');
   wbdf = get(disph,'WindowbuttonDownFcn');
   set(disph,'WindowbuttonDownFcn','');
   axes(ax);
   initpt = get(ax,'currentpoint');
   initpt = initpt(1,[1,2]);
   pplane5('solution',initpt,disph);
   set(disph,'WindowbuttonDownFcn',wbdf);
   %  set([ch;mh],'enable','on');
   set(ch,'enable','on');
   notice = dud.notice;
   if notice
      nstr = get(notice,'string');
      nstr(1:4) = nstr(2:5);
      nstr{5,1} = 'Ready.';
      set(notice,'string',nstr)
   end
   
elseif strcmp(action,'several')
   
   % 'several' allows the user to pick several initial points at once.
   % This is not needed in X-windows, but it is on the Macintosh.
   
   disph = gcf;
   ch = findobj('type','uicontrol','enable','on');
   %  mh = findobj('type','uimenu','enable','on');
   set(ch,'enable','inactive');
   %  set(mh,'enable','off')
   wbdf = get(disph,'WindowbuttonDownFcn');
   set(disph,'WindowbuttonDownFcn','');
   dud = get(disph,'user');
   notice = dud.notice;
   if notice
      nstr = get(notice,'string');
      nstr(1:4) = nstr(2:5);
      nstr{5,1} = 'Pick initial points with the mouse. Enter "Return" when finished.';
      set(notice,'string',nstr)
   end
   [X,Y]=ginput;
   NN = length(X);
   for k = 1:NN
      initpt = [X(k),Y(k)];
      pplane5('solution',initpt,disph);
   end
   if notice
      nstr = get(notice,'string');
      nstr(1:4) = nstr(2:5);
      nstr{5,1} = 'Ready.';
      set(notice,'string',nstr);
   end
   set(ch,'enable','on');
   %  set([ch;mh],'enable','on');
   set(disph,'WindowbuttonDownFcn',wbdf);
   
elseif strcmp(action,'test case')
   
   tic  
   ppdisp = gcf;
   for k = 1:10
      initpt = k*[-.2 .2];
      pplane5('solution',initpt,ppdisp);
      initpt = k*[-.2 -.2];
      pplane5('solution',initpt,ppdisp);
   end  
   pplane5('solution',[1,-1],ppdisp);
   dud = get(gcf,'user');
   nstr = get(dud.notice,'string');
   nstr(1:4) = nstr(2:5);
   nstr{5,1} = 'Ready.';
   set(dud.notice,'string',nstr)
   toc
   
elseif strcmp(action,'solution')
   
   % 'solution' effects the computation and (erasemode == none) plotting of
   % solutions.  It also stores the data as appropriate.
   
   disph = input2;
   dud = get(disph,'user');
   tcol = dud.color.temp;
   pcol = dud.color.orb;
   notice = dud.notice;
   initpt = input1(:);
   dfcn = dud.function;
   ppdispa = dud.axes;
   settings = dud.settings;
   ptstr = [' (',num2str(initpt(1),2), ', ', num2str(initpt(2),2), ')'];
   refine = settings.refine;
   ssize = settings.stepsize;
   tol = settings.tol;
   ud = get(dud.axes,'user');
   %  rtol = tol;
   atol = tol*ud.DY*1e-4';
   
   
   if length(initpt)  == 2
      AA = -1e6;
      BB = 1e6; 
      switch dud.dir
         %      case 0
         %	intplus = [initpt(1), BB];
         %	intminus = [initpt(1), AA]; 
         %      case -1
         %	intplus = [initpt(1), initpt(1)];
         %	intminus = [initpt(1), AA]; 
         %      case 1
         %	intplus = [initpt(1), BB];
         %	intminus = [initpt(1), initpt(1)];
      case 0
         intplus = [0, BB];
         intminus = [0, AA]; 
      case -1
         intplus = [0, 0];
         intminus = [0, AA]; 
      case 1
         intplus = [0, BB];
         intminus = [0, 0];
      end
      
   else
      intplus = [initpt(3),initpt(5)];
      intminus = [initpt(3),initpt(4)];
      initpt = initpt([1:2]);
      
   end	    
   stopbutt = findobj(disph,'tag','stop');
   
   set(stopbutt,'vis','on','enable','on');
   
   options = odeset('OutputFcn','ppn5out','Refine',refine,...
       'RelTol',tol,'Abstol',atol);

   cflag = 0;
   if intplus(2)>intplus(1)
      cflag = cflag + 1;
      if notice
         nstr = get(notice,'string');
         nstr(1:4) = nstr(2:5);
         nstr{5} = ['The forward orbit from',ptstr];
         set(notice,'string',nstr);
      end
      drawnow				
      
      switch settings.solver
	case 'Dormand Prince'	
	  [tp,xp] = ppdp45(dfcn,intplus,initpt,disph);
	  
	case 'ode45'	
	  [tp,xp] = ode45(dfcn,intplus,initpt,options);
	  
	case 'ode23'	
	  [tp,xp] = ode23(dfcn,intplus,initpt,options);
	  
	case 'ode113'	
	  [tp,xp] = ode113(dfcn,intplus,initpt,options);
	  
	case 'ode15s'	
	  [tp,xp] = ode15s(dfcn,intplus,initpt,options);
	  
	case 'ode23s'	
	  [tp,xp] = ode23s(dfcn,intplus,initpt,options);
	  
	case 'Runge-Kutta 4'	
	  [tp,xp] = pprk4(dfcn,intplus,initpt,disph);
	  
	otherwise
	  error('Undefined solver.');
      end	  

      aud = get(ppdispa,'user');
      hnew1 = aud.line;
   end		
   
   if intminus(2) < intminus(1)
      cflag = cflag + 2;
      if notice	
         nstr = get(notice,'string');
         nstr(1:4) = nstr(2:5);
         nstr{5} = ['The backward orbit from',ptstr];
         set(notice,'string',nstr);
      end
      drawnow
      switch settings.solver
	case 'Dormand Prince'
	  [tm,xm] = ppdp45(dfcn,intminus,initpt,disph);
	  
	case 'ode45'	
	  [tm,xm] = ode45(dfcn,intminus,initpt,options);
	  
	case 'ode23'	
	  [tm,xm] = ode23(dfcn,intminus,initpt,options);
	  
	case 'ode113'	
	  [tm,xm] = ode113(dfcn,intminus,initpt,options);
	  
	case 'ode15s'	
	  [tm,xm] = ode15s(dfcn,intminus,initpt,options);
	  
	case 'ode23s'	
	  [tm,xm] = ode23s(dfcn,intminus,initpt,options);
	  
	case 'Runge-Kutta 4'	
	  [tm,xm] = pprk4(dfcn,intminus,initpt,disph);
	  
	otherwise
	  error('Undefined solver.');
      end

      aud = get(ppdispa,'user');
      hnew2 = aud.line;
      
      set(stopbutt,'vis','off');
   end	  % if intminus(2) < intminus(1)	
   
   
   
   % Store the trajectory.
   
   switch cflag
   case 1 % positive only
      set(hnew1,'xdata',xp(:,1),'ydata',xp(:,2),'zdata',tp,'color',pcol);
      set(hnew1,'erase','normal');
      dud.solhand = [dud.solhand;hnew1];
   case 2 % negative only
      x = flipud(xm);
      t = flipud(tm);
      set(hnew2,'xdata',x(:,1),'ydata',x(:,2),'zdata',t,'color',pcol);
      set(hnew2,'erase','normal');
      dud.solhand = [dud.solhand;hnew2];
      
   case 3 % both directions
      x = flipud(xm);
      t = flipud(tm);
      x=[x;xp];
      t=[t;tp];
      delete(hnew1);
      set(hnew2,'xdata',x(:,1),'ydata',x(:,2),'zdata',t,'color',pcol);
      set(hnew2,'erase','normal');
      dud.solhand = [dud.solhand;hnew2];
   end	 % switch cflag
   set(disph,'user',dud);
   
elseif strcmp(action,'kcompute')
   
   % 'kcompute' is the call back for the Compute 
   % button on the PPLANE5 Keyboard figure.
   
   compute = 1;
   kh = get(gcf,'user');
   ppdisp = findobj('name','PPLANE5 Display');
   if (isempty(ppdisp))
      pplane5('confused');
   end
   dud = get(ppdisp,'user');
   ppdispa = dud.axes;
   aud = get(ppdispa,'user');
   ch = findobj('type','uicontrol','enable','on');
   %  mh = findobj('type','uimenu','enable','on');
   set(ch,'enable','inactive');
   %  set(mh,'enable','off')
   set(ppdisp,'WindowbuttonDownFcn','');
   x = str2num(get(kh.xval,'string'));
   y = str2num(get(kh.yval,'string'));
   if get(kh.spec,'value')
      tzero = str2num(get(kh.tval,'string'));
      t0 = str2num(get(kh.t0,'string'));
      tf = str2num(get(kh.tf,'string'));
      initpt = [x,y,tzero,t0,tf];
      if (length(initpt) ~= 5)
         warndlg({'Values must be entered for all of the entries.'},'Illegal input');
         compute = 0;
      elseif tf <= t0
         warndlg({'The final time of the computation interval';...
               'must be smaller than the initial time.'},'Illegal input');
         compute = 0;
      elseif (tzero < t0) | (tzero > tf)
         warndlg('The initial time must be in the computation interval.','Illegal input');
         compute = 0;
      end
      aud.gstop = 0;
      set(ppdispa,'user',aud);
   else
      initpt = [x,y];
      if (length(initpt) ~= 2)
         warndlg({'Values must be entered for all of the entries.'},'Illegal input');
         compute = 0;
      end
   end	% if get(kh.spec,'value')
   
   if compute
      pplane5('solution',initpt,ppdisp);
   end
   nstr = get(dud.notice,'string');
   nstr(1:4) = nstr(2:5);
   nstr{5} = 'Ready.';
   set(dud.notice,'string',nstr)
   %  set([ch;mh],'enable','on');
   set(ch,'enable','on');
   set(ppdisp,'WindowbuttonDownFcn','pplane5(''down'')');
   aud.gstop = 1;
   set(ppdispa,'user',aud);
   
   
elseif strcmp(action,'kbd')
   
   % 'kbd' is the callback for the Keyboard Input menu selection.  It 
   % sets up the PPLANE5 Keyboard figure which allows accurate input of 
   % initial values using the keyboard.
   
   ppkbd = findobj('name','PPLANE5 Keyboard input');
   ppdisp = gcf;
   dud = get(ppdisp,'user');
   Xname = dud.syst.xvar;
   Yname = dud.syst.yvar;
   xnstr = ['The initial value of ',Xname,' = '];
   ynstr = ['The initial value of ',Yname,' = '];
   tnstr = 'The initial value of t = ';
   nudge = 3; 
   if  isempty(ppkbd)
      
      ppkbd = figure('name','PPLANE5 Keyboard input',...
         'vis','off',...
         'numb','off','tag','PPLANE5');
      
      pplane5('figdefault',ppkbd);
      fmenu = findobj(gcf,'label','File');  % uimenu('label','File');
      clmenu = uimenu(fmenu,'label','Close',...
         'call','close','accelerator','W');
      compmenu = uimenu(fmenu,'label','Compute',...
         'call','pplane5(''kcompute'')','accelerator','C');
      qmenu = uimenu(fmenu,'label','Quit PPLANE5',...
         'call','pplane5(''quit'')','separator','on');
      
      kbd.fr1 = uicontrol('style','frame');
      
      kbd.fr2 = uicontrol('style','frame');
      
      kbd.inst = uicontrol('style','text','horiz','center',...
         'string','Enter the initial conditions:');
      
      kbd.xname = uicontrol('style','text',...
         'horiz','right','string',xnstr);
      
      kbd.xval = uicontrol('style','edit',...
         'string','','call','');
      
      kbd.yname = uicontrol('style','text',...
         'horiz','right',...
         'string',ynstr);
      
      kbd.yval = uicontrol('style','edit',...
         'string','');
      
      kbd.tname = uicontrol('style','text',...
         'horiz','right',...
         'string',tnstr);
      
      kbd.tval = uicontrol('style','edit');
      
      kbd.t0 = uicontrol('style','edit');
      
      kbd.tf = uicontrol('style','edit');
      
      kbd.t = uicontrol('style','text','string','<= t<= ');
      
      kbd.spec = uicontrol('style','check','horiz','center',...
         'string','Specify a computation interval.');
      
      kbd.comp = uicontrol('style','push',...		
         'string','Compute','call','pplane5(''kcompute'')');
      
      kbd.close = uicontrol('style','push',...
         'string','Close','call','set(gcf,''vis'',''off'')');
      
      left = 5; varl = 70; buttw = 60;
      xex = get(kbd.xname,'extent');
      ht = xex(4)+nudge;
      yex = get(kbd.yname,'extent');
      nwdth = max(xex(3),yex(3)) + nudge;
      varl = varl*ht/19;
      fr1bot = 2*left + ht;
      fr1ht = 3*nudge + 2*ht;
      fr2bot = fr1bot + fr1ht;
      frw = 2*nudge + nwdth + varl;
      intb = fr1bot + nudge;
      specb = intb + ht + nudge;
      intext = get(kbd.t,'extent');
      tw = intext(3);
      margin = (frw -tw - 2*varl)/2;
      t0l = left+margin;
      tl = t0l + varl;
      tfl = tl + tw;
      tnb = fr2bot + nudge;
      ynb = tnb + ht;
      xnb = ynb + ht;
      nl = left+nudge;
      xnwind = [nl,xnb,nwdth,ht];
      ynwind = [nl,ynb,nwdth,ht];
      tnwind = [nl,tnb,nwdth,ht];
      vl = nl + nwdth;
      xvwind = [vl,xnb,varl,ht];
      yvwind = [vl,ynb,varl,ht];
      tvwind = [vl,tnb,varl,ht];
      instb = xnb + ht + nudge;
      instw = nwdth + varl;
      instwind = [nl,instb,instw,ht];
      fr2ht = 3*nudge + 4*ht;
      frw = 2*nudge + nwdth + varl;
      fr2wind = [left,fr2bot,frw,fr2ht];
      fr1wind = [left,fr1bot,frw,fr1ht];
      specw = frw - 2*nudge;
      specwind = [nl,specb,specw,ht];
      t0wind = [t0l,intb,varl,ht];
      twind = [tl,intb,tw,ht];
      tfwind = [tfl,intb,varl,ht];
      figw = 2*left + frw;
      fight = 3*left + ht + fr1ht + fr2ht;
      figwind = [30,300,figw,fight];
      buttw = frw/3;
      sep = (figw - 2*buttw)/3;
      closel = sep;
      compl = 2*sep+buttw;
      clwind = [closel,left,buttw,ht];
      compwind = [compl,left,buttw,ht];
      set(ppkbd,'pos',figwind);
      set(kbd.fr1,'pos',fr1wind);
      set(kbd.fr2,'pos',fr2wind);
      set(kbd.inst,'pos',instwind);
      set(kbd.xname,'pos',xnwind);
      set(kbd.yname,'pos',ynwind);
      set(kbd.xval,'pos',xvwind);
      set(kbd.yval,'pos',yvwind);
      set(kbd.comp,'pos',compwind);
      set(kbd.close,'pos',clwind);
      set(kbd.tname,'pos',tnwind);
      set(kbd.tval,'pos',tvwind);
      set(kbd.spec,'pos',specwind);
      set(kbd.t0,'pos',t0wind);
      set(kbd.t,'pos',twind);
      set(kbd.tf,'pos',tfwind);
      speccall = [
         'ud = get(gcf,''user'');',...
            'if get(gcbo,''value''),',...
            '  set([ud.t0,ud.t,ud.tf,ud.tname,ud.tval],''vis'',''on'');',...
            'else,',...
            '  set([ud.t0,ud.t,ud.tf,ud.tname,ud.tval],''vis'',''off'');',...
            'end'];
      
      set(kbd.spec,'call',speccall);
      set(ppkbd,'resize','on');
      set(findobj(ppkbd,'type','uicontrol'),'units','normal');
      
   else
      set(findobj(ppkbd,'type','uicontrol'),'units','pix');
      kbd = get(ppkbd,'user');
      set(kbd.xname,'string',xnstr);
      set(kbd.yname,'string',ynstr);
      set(kbd.xval,'string','');
      set(kbd.yval,'string','');
      xpos = get(kbd.xname,'pos');
      xex = get(kbd.xname,'extent');
      yex = get(kbd.yname,'extent');
      nwdth = max(xex(3),yex(3)) + nudge;
      delta = nwdth - xpos(3);
      if delta ~= 0
         xpos(3) = nwdth;
         set(kbd.xname,'pos',xpos);
         set(kbd.yname,'pos',get(kbd.yname,'pos')+[0,0,delta,0]);
         set(kbd.tname,'pos',get(kbd.tname,'pos')+[0,0,delta,0]);
         set(kbd.xval,'pos',get(kbd.xval,'pos')+[delta,0,0,0]);
         set(kbd.yval,'pos',get(kbd.yval,'pos')+[delta,0,0,0]);
         set(kbd.tval,'pos',get(kbd.tval,'pos')+[delta,0,0,0]);
         set(kbd.fr1,'pos',get(kbd.fr1,'pos')+[0,0,delta,0]);
         set(kbd.fr2,'pos',get(kbd.fr2,'pos')+[0,0,delta,0]);
         set(kbd.spec,'pos',get(kbd.spec,'pos')+[0,0,delta,0]);
         set(kbd.inst,'pos',get(kbd.inst,'pos')+[0,0,delta,0]);
         set(ppkbd,'pos',get(ppkbd,'pos')+[0,0,delta,0]);
         set(kbd.t0,'pos',get(kbd.t0,'pos')+[delta/2,0,0,0]);
         set(kbd.t,'pos',get(kbd.t,'pos')+[delta/2,0,0,0]);
         set(kbd.tf,'pos',get(kbd.tf,'pos')+[delta/2,0,0,0]);
         set(kbd.comp,'pos',get(kbd.comp,'pos')+[2*delta/3,0,0,0]);
         set(kbd.close,'pos',get(kbd.close,'pos')+[delta/3,0,0,0]);
      end % delta ~= 0
   end % if  isempty(ppkbd) & else
   set([kbd.tval,kbd.t0],'string','0');
   set(kbd.spec,'value',0);
   set(ppkbd,'user',kbd,'vis','on');
   set([kbd.t0,kbd.t,kbd.tf,kbd.tname,kbd.tval],'vis','off')
   set(findobj(ppkbd,'type','uicontrol'),'units','normal');
   figure(ppkbd)
   
elseif strcmp(action,'eqpt')
   
   % Find and classify equilibrium points.
   
   ppdisp = findobj('name','PPLANE5 Display');
   dud = get(ppdisp,'user');
   col = dud.color.eqpt;
   dfcn = dud.function;
   dbutt = dud.butt;
   menu = dud.menu;
   ppdispa =dud.axes;
   Dx = get(ppdispa,'xlim');
   Dy = get(ppdispa,'ylim');
   DY = [Dx(2)-Dx(1);Dy(2)-Dy(1)];
   
   epsilon = 1e-4;
   nstr = get(dud.notice,'string');
   nstr(1:4) = nstr(2:5);
   nstr{5} = 'Choose an approximation with the mouse.';
   set(dud.notice,'string',nstr);
   
   z0 = ginput(1);
   
   z = pplane5('newton',z0,dfcn);
   
   flag = z(:,4);
   A = z(:,2:3);
   B = real(A);
   k = find(abs(B)<1e-8);
   B(k) = zeros(size(k));
   C = imag(A);
   k = find(abs(C)<1e-8);
   C(k) = zeros(size(k));
   A = B+C*sqrt(-1);
   z = z(:,1);
   
   if (~flag | norm((z-z0')./DY) > 1/5)
      nstr(1:4) = nstr(2:5);
      nstr{5} = ['There is not an equilibrium point near (',...
            num2str(z0(1)),', ', num2str(z0(2)), ').  Ready'];
      set(dud.notice,'string',nstr);
      return
   end
   zero =find(abs(z) < epsilon);
   z(zero) = zeros(size(zero));
   D=det(A); T=trace(A);
   if (D<-epsilon)
      string1 = ['There is a saddle point at (',...
            num2str(z(1)),', ', num2str(z(2)), ').'];
      string2 = '';
      type = 1;
   elseif(D>epsilon)
      if (T*T>4*D+epsilon)
         if (T<0)
            string1 = ['There is a nodal sink at (',...
                  num2str(z(1)),', ', num2str(z(2)), ').'];
            string2 = '';
            type = 2;
         else
            string1 = ['There is a nodal source at (',...
                  num2str(z(1)),', ', num2str(z(2)), ').'];
            string2 = '';
            type = 3;
         end
      elseif(T*T<4*D-epsilon)
         if(T<-epsilon)
            string1 = ['There is a spiral sink at (',...
                  num2str(z(1)),', ', num2str(z(2)), ').'];
            string2 = '';
            type = 4;
         elseif(T>epsilon)
            string1 = ['There is a spiral source at (',...
                  num2str(z(1)),', ', num2str(z(2)), ').'];
            string2 = '';
            type = 5;
         else
            string1 = ['There is a spiral equilibrium point at (',...
                  num2str(z(1)),', ', num2str(z(2)), ').'];
            string2 = 'Its specific type has not been determined.';
            type = 6;
         end
      else
         if (T>epsilon)
            string1 = ['There is a source at (',...
                  num2str(z(1)),', ', num2str(z(2)), ').'];
            string2 = 'Its specific type has not been determined.';
            type = 7;
         elseif (T<-epsilon)
            string1 = ['There is a sink at (',...
                  num2str(z(1)),', ', num2str(z(2)), ').'];
            string2 = 'Its specific type has not been determined.';
            type = 8;
         else
            string1 = ['There is an equilibrium point at (',...
                  num2str(z(1)),', ', num2str(z(2)), ').'];
            string2 = 'Its specific type has not been determined.';
            type = 9;
         end
      end
   else
      string1 = ['There is an equilibrium point at (',...
            num2str(z(1)),', ', num2str(z(2)), ')'];
      string2 = 'Its specific type has not been determined.';
      type = 9;
   end
   
   [V,D] = eig(A);			
   EqPtList = dud.eqpts;
   infostr{1,1} = string1;
   infostr{2,1} = string2;
   infostr{3,1} = ' ';
   infostr{4,1} = 'The Jacobian is:';
   infostr{5,1} = ['      ',num2str(A(1,1)),'   ',num2str(A(1,2))];
   infostr{6,1} = ['      ',num2str(A(2,1)),'   ',num2str(A(2,2))];
   infostr{7,1} = '';
   infostr{8,1} = 'The eigenvalues and eigenvectors are:';
   infostr{9,1} = ['     ',num2str(D(1,1)),'    (',num2str(V(1,1)),', ',num2str(V(2,1)),')'];
   infostr{10,1} = ['     ',num2str(D(2,2)),'    (',num2str(V(1,2)),', ',num2str(V(2,2)),')'];
   
   k=1; l=size(EqPtList,1);
   while ((k <= l))
      if (norm((EqPtList(k,1:2)-z')./DY') <= 1e-3)
         break;
      end
      k = k+1;
   end
   if (k > l)
      EqPtList = [EqPtList;z',type];
      dud.eqpts = EqPtList;
      newh = plot(z(1),z(2),...
         'markerfacecolor',col,...
         'markeredgecolor',col,...
         'marker','o','Erasemode','none');
      dud.solhand = [dud.solhand;newh(:)];
      set(newh,'Erasemode','normal');
      drawnow
   end
   ppeqpt = findobj('name','PPLANE5 Equilibrium point data');
   
   if (isempty(ppeqpt))
      ppeqpt = figure('vis','off','resize','on',...
         'name','PPLANE5 Equilibrium point data',...
         'numb','off','tag','PPLANE5');
      pplane5('figdefault',ppeqpt);
      fmenu = findobj(ppeqpt,'label','File');   % uimenu('label','File');
      clmenu = uimenu(fmenu,'label','Close',...
         'call','close',...
         'accelerator','W');
      qmenu = uimenu(fmenu,'label','Quit PPLANE5',...
         'call','pplane5(''quit'')',...
         'separator','on');
      ud.frame = uicontrol('style','frame');
      ud.eptext = uicontrol('style','text','string',infostr,'hor','left');
      ud.goaway = uicontrol('style','push','string','Go away',...
         'call','close');
      ud.display = uicontrol('style','push',...
         'string','Display the linearization',...
         'call','pplane5(''linear'')');
   else
      figure(ppeqpt);
      ud = get(ppeqpt,'user');
      set(ppeqpt,'vis','off');
      set(findobj(ppeqpt,'type','uicontrol'),'units','pixels');
      set(ud.eptext,'string',infostr);
   end
   ud.jac = A;
   ud.vectors = V;
   ud.type = type;
   ud.system = dud.syst;
   ud.settings = dud.settings;
   ud.color = dud.color;
   set(ppeqpt,'user',ud);
   left = 5; nudge = 3;
   ext = get(ud.eptext,'extent');
   n = size(infostr,1);
   ht = ext(4)/n+2*nudge;
   frbot = 3*left + 2*ht;
   txtbot = frbot + nudge;
   txtl = left + nudge;
   twind = [txtl,txtbot,ext(3)+nudge,ext(4)];
   frw = ext(3) + 3*nudge;
   frh = ext(4) + 2*nudge;
   frwind = [left,frbot,frw,frh];
   figw = frw + 2*left;
   figh = frh + 4*left + 2*ht;
   uni = get(0,'units');
   set(0,'units','pixels');
   ss = get(0,'screensize');
   set(0,'units',uni);
   sh = ss(4);
   figbot = sh - figh - 350;
   figwind = [30,figbot,figw,figh];
   buttw = frw - 4;
   buttl = left + 2;
   closewind = [buttl,left,buttw,ht];
   dispb = 2*left + ht;
   dispwind = [buttl,dispb,buttw,ht];
   set(ppeqpt,'pos',figwind);
   set(ud.frame,'pos',frwind);
   set(ud.eptext,'pos',twind);
   set(ud.goaway,'pos',closewind);
   set(ud.display,'pos',dispwind);
   set(findobj(ppeqpt,'type','uicontrol'),'units','normal');
   set(ppeqpt,'vis','on');
   set(get(ppeqpt,'child'),'vis','on');
   nstr(1:4) = nstr(2:5);
   nstr{5} = 'Ready.';
   set(dud.notice,'string',nstr);
   set(ppdisp,'user',dud);
   
elseif strcmp(action,'stunst')
   
   ppdisp = findobj('name','PPLANE5 Display');
   dud = get(ppdisp,'user');
   ecol = dud.color.eqpt;
   scol = dud.color.sep;
   ppdispa = dud.axes;
   aud = get(ppdispa,'user');
   DY = aud.DY;
   settings = dud.settings;
   dfcn = dud.function;
   EqPtList = dud.eqpts;
   Stop = norm(DY)*1e-4;
   
   %  Plot the stable and unstable orbits at a saddle point.
   
   nstr = get(dud.notice,'string');
   nstr(1:4) = nstr(2:5);
   nstr{5} = 'Choose a saddle point with the mouse.';
   set(dud.notice,'string',nstr);
   z0 = ginput(1);
   z = zeros(2,1);
   z = pplane5('newton',z0,dfcn);
   
   flag = z(:,4);
   A = z(:,2:3);
   k = find(abs(A)<1e-8);
   A(k) = zeros(size(k));
   z = z(:,1);
   
   if (~flag | norm((z-z0')./DY)> 1/5)
      nstr(1:4) = nstr(2:5);
      nstr{5} = ['There is not an equilibrium point near (',...
            num2str(z0(1),2),', ', num2str(z0(2),2), ').'];
      set(dud.notice,'string',nstr);
      return
   end
   zero = find(abs(z) < 1e-4);
   z(zero) = zeros(size(zero));
   D=det(A);
   if (D>=0)
      nstr(1:4) = nstr(2:5);
      nstr{5} = ['The equilibrium point at (',...
            num2str(z(1),2),', ', num2str(z(2),2),...
            ') is not a saddle point.'];
      set(dud.notice,'string',nstr);
      return
   end
   if isempty(EqPtList)
      EqPtList = [z',1];
      dud.eqpts = EqPtList;
   else
      k=1; l=size(EqPtList,1);
      while ((k <= l))
         if (norm(EqPtList(k,1:2)-z') <= Stop)
            break;
         end
         k = k+1;
      end
      if (k > l)
         EqPtList = [EqPtList;z',1];
         dud.eqpts = EqPtList;
      end
   end
   nstr(1:4) = nstr(2:5);
   nstr{5} = ['Plotting from the saddle point at (',...
         num2str(z(1),2),', ', num2str(z(2),2),').'];
   set(dud.notice,'string',nstr);
   [V,L] = eig(A);
   [L,I] = sort(diag(L));
   
   magn = settings.magn;
   refine = settings.refine;
   ssize = settings.stepsize;
   solver = settings.solver;
   tol = settings.tol;
   offset = norm(DY)/1000;	
   lm=abs(L);
   lm=lm/max(lm);
   lm = (max([lm,0.2*ones(size(lm))]'))';
   offs = offset./lm;
   ud = get(dud.axes,'user');
   atol = tol*ud.DY*1e-4';
   options = odeset('OutputFcn','ppn5out','Refine',refine,...
       'RelTol',tol,'Abstol',atol);
   stopbutt = findobj('tag','stop');
   set(stopbutt,'vis','on','enable','on');
   VV=V(:,I(1));
   hnew = zeros(5,1);
   w = z + offs(1)*VV;
   nstr = get(dud.notice,'string');
   nstr(1:4) = nstr(2:5);
   nstr{5} = 'The first stable trajectory';
   set(dud.notice,'string',nstr);
   int = [0,-1e6];

   switch settings.solver
     case 'Dormand Prince'	
       [tp,Xp] = ppdp45(dfcn,int,w,ppdisp);
       
     case 'ode45'	
       [tp,Xp] = ode45(dfcn,int,w,options);
       
     case 'ode23'	
       [tp,Xp] = ode23(dfcn,int,w,options);
       
     case 'ode113'	
       [tp,Xp] = ode113(dfcn,int,w,options);
       
     case 'ode15s'	
       [tp,Xp] = ode15s(dfcn,int,w,options);
       
     case 'ode23s'	
       [tp,Xp] = ode23s(dfcn,int,w,options);
       
     case 'Runge-Kutta 4'	
       [tp,Xp] = pprk4(dfcn,int,w,ppdisp);
       
     otherwise
       error('Undefined solver.');
   end	  
%   [tp,Xp] = ppdp45(dfcn,[0,-1e6],w,ppdisp);
   set(stopbutt,'enable','off');
   aud = get(ppdispa,'user');
   hnew(1) = aud.line;
   X1=[z';Xp];
   x1 = [X1,[NaN;tp]];
   set(hnew(1),'xdata',x1(:,1),'ydata',x1(:,2),'zdata',x1(:,3));
   
   w = z - offs(1)*VV;
   nstr = get(dud.notice,'string');
   nstr(1:4) = nstr(2:5);
   nstr{5} = 'The second stable trajectory';
   set(dud.notice,'string',nstr);
   set(stopbutt,'enable','on');

   switch settings.solver
     case 'Dormand Prince'	
       [tp,Xp] = ppdp45(dfcn,int,w,ppdisp);
       
     case 'ode45'	
       [tp,Xp] = ode45(dfcn,int,w,options);
       
     case 'ode23'	
       [tp,Xp] = ode23(dfcn,int,w,options);
       
     case 'ode113'	
       [tp,Xp] = ode113(dfcn,int,w,options);
       
     case 'ode15s'	
       [tp,Xp] = ode15s(dfcn,int,w,options);
       
     case 'ode23s'	
       [tp,Xp] = ode23s(dfcn,int,w,options);
       
     case 'Runge-Kutta 4'	
       [tp,Xp] = pprk4(dfcn,int,w,ppdisp);
       
     otherwise
       error('Undefined solver.');
   end	  
%   [tp,Xp] = ppdp45(dfcn,[0,-1e6],w,ppdisp);
   set(stopbutt,'enable','off');
   
   aud = get(ppdispa,'user');
   hnew(2) = aud.line;
   
   X2=[z';Xp];
   x2 = [X2,[NaN;tp]];
   set(hnew(2),'xdata',x2(:,1),'ydata',x2(:,2),'zdata',x2(:,3));
   
   VV=V(:,I(2));
   w = z + offs(2)*VV;
   nstr = get(dud.notice,'string');
   nstr(1:4) = nstr(2:5);
   nstr{5} = 'The first unstable trajectory';
   set(dud.notice,'string',nstr);
   set(stopbutt,'enable','on');
   int = [0,1e6];

   switch settings.solver
     case 'Dormand Prince'	
       [tp,Xp] = ppdp45(dfcn,int,w,ppdisp);
       
     case 'ode45'	
       [tp,Xp] = ode45(dfcn,int,w,options);
       
     case 'ode23'	
       [tp,Xp] = ode23(dfcn,int,w,options);
       
     case 'ode113'	
       [tp,Xp] = ode113(dfcn,int,w,options);
       
     case 'ode15s'	
       [tp,Xp] = ode15s(dfcn,int,w,options);
       
     case 'ode23s'	
       [tp,Xp] = ode23s(dfcn,int,w,options);
       
     case 'Runge-Kutta 4'	
       [tp,Xp] = pprk4(dfcn,int,w,ppdisp);
       
     otherwise
       error('Undefined solver.');
   end	  
%   [tp,Xp] = ppdp45(dfcn,[0,1e6],w,ppdisp);
   set(stopbutt,'enable','off');
   
   aud = get(ppdispa,'user');
   hnew(3) = aud.line;
   X3=[z';Xp];
   x3 = [X3,[NaN;tp]];
   set(hnew(3),'xdata',x3(:,1),'ydata',x3(:,2),'zdata',x3(:,3));
   
   w = z - offs(2)*VV;
   nstr = get(dud.notice,'string');
   nstr(1:4) = nstr(2:5);
   nstr{5} = 'The second unstable trajectory';
   set(dud.notice,'string',nstr);
   set(stopbutt,'enable','on');
   
   switch settings.solver
     case 'Dormand Prince'	
       [tp,Xp] = ppdp45(dfcn,int,w,ppdisp);
       
     case 'ode45'	
       [tp,Xp] = ode45(dfcn,int,w,options);
       
     case 'ode23'	
       [tp,Xp] = ode23(dfcn,int,w,options);
       
     case 'ode113'	
       [tp,Xp] = ode113(dfcn,int,w,options);
       
     case 'ode15s'	
       [tp,Xp] = ode15s(dfcn,int,w,options);
       
     case 'ode23s'	
       [tp,Xp] = ode23s(dfcn,int,w,options);
       
     case 'Runge-Kutta 4'	
       [tp,Xp] = pprk4(dfcn,int,w,ppdisp);
       
     otherwise
       error('Undefined solver.');
   end	  
%   [tp,Xp] = ppdp45(dfcn,[0,1e6],w,ppdisp);
   set(stopbutt,'vis','off');
   
   aud = get(ppdispa,'user');
   hnew(4) = aud.line;
   X4=[z';Xp];
   x4 = [X4,[NaN;tp]];
   set(hnew(4),'xdata',x4(:,1),'ydata',x4(:,2),'zdata',x4(:,3));
   
   hnew(5) = plot(z(1),z(2),...
      'markerfacecolor',ecol,...
      'markeredgecolor',ecol,...
      'marker','o','Erasemode','none');
   %  drawnow
   dud.solhand = [dud.solhand;hnew(:)];
   set(hnew(1:4),'color',scol);
   set(hnew,'erase','normal');
   nstr = get(dud.notice,'string');
   nstr(1:4) = nstr(2:5);
   nstr{5} = 'Ready.';
   set(dud.notice,'string',nstr);
   set(gcf,'user',dud);
   
elseif strcmp(action,'zoomin')
   
   % 'zoomin' is the callback for the Zoomin menu item.  It allows the
   % user to pick a new display rectangle.
   
   set(gcf,'WindowButtonDownFcn','pplane5(''zoom'')',...
      'WindowButtonUpFcn','1;','inter','on');
   set(gca,'inter','on');
   dud = get(gcf,'user');
   nstr = get(dud.notice,'string');
   nstr(1:4) = nstr(2:5);
   nstr{5} = ['Pick a new display rectangle by clicking and ',...
         'dragging the mouse, or by clicking on a point.'];
   set(dud.notice,'string',nstr);
   
elseif strcmp(action,'zoom')
   
   disph = gcf;
   dud = get(disph,'user');
   axh = dud.axes;
   aud = get(axh,'user');
   DY = aud.DY;
   w = dud.syst.wind;
   q1 = get(disph,'currentpoint');
   p1 = get(axh,'currentpoint');
   p1 = p1(1,1:2);
   rbbox([q1 0 0],q1);
   p2 = get(axh,'currentpoint');
   p2 = p2(1,1:2);
   if all(abs(p1'-p2')>0.01*DY)
      a = [p1;p2];
      a = [min(a);max(a)];
      DY = (a(2,:) - a(1,:))';  
   else
      DY = DY/4;
      a(1) = max(w(1),p1(1)-DY(1));
      a(2) = min(w(2),p1(1)+DY(1));
      a(3) = max(w(3),p1(2)-DY(2));
      a(4) = min(w(4),p1(2)+DY(2));
      DY(1) = a(2) - a(1);
      DY(2) = a(4) - a(3);
   end    
   WINvect = a(:)';
   dud.syst.wind = WINvect;
   aud.DY = DY;
   dwind = [WINvect(1); WINvect(3); -WINvect(2); -WINvect(4)];
   aud.cwind = dwind - dud.settings.magn*[aud.DY;aud.DY];
   set(axh,'user',aud);	
   set(disph,'user',dud);
   set(disph,'WindowButtonDownFcn','pplane5(''down'')',...
      'WindowButtonUpFcn','');
   pplane5('dirfield',disph);
   ppset = findobj('name','PPLANE5 Setup');
   if isempty(ppset)
      pplane5('confused');
   else
      sud = get(ppset,'user');
      sud.c.wind = WINvect;
      sud.o.wind = WINvect;
      set(sud.h.wind(1),'string',num2str(WINvect(1)));
      set(sud.h.wind(2),'string',num2str(WINvect(2)));
      set(sud.h.wind(3),'string',num2str(WINvect(3)));
      set(sud.h.wind(4),'string',n