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

% PPLANE  is an interactive tool for studying planar autonomous systems of
%	differential equations.  When PPLANE is executed, a PPLANE 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 PPLANE
%	Display window is opened, and a field of the type requested is 
%	displayed for the system.  When the mouse button is depressed in 
%	the PPLANE Display window, the solution to the system with that i
%	nitial condition is calculated and plotted.
%
%	Other options are available in the Options menu.  These are
%	fairly self explanatory.  The Settings option allows the user
%	to change several parameters affecting the speed and accuracy of 
%	the computation.

%  Copywright (c)  John C. Polking, Rice University
%  Last Modified: September 24, 1995


if nargin <1
  action ='initialize';
end


if strcmp(action,'initialize')

  % First we make sure that there is no other copy of PPLANE
  % running, since this causes problems.	

  figs = get(0,'children');
  ppdisp = 0;
  for (figno = 1:length(figs))
    nn = get(figs(figno),'name');
    if (length(nn) > 3)
      if (findstr(nn,'PPLANE')~=[])
	error('Only one copy of PPLANE can be open at one time.')
      end
    end
  end
	
  % Set up for the menu of systems.

  default1 = str2mat(...
      'x','y',...
      ' 2*x - y + 3*(x^2-y^2) + 2*x*y',...
      ' x - 3*y - 3*(x^2-y^2) + 3*x*y',...
      ' ',' ',' ',' ');
  default2 = str2mat(...
      ' ',' ',' ','  ',...
      '2','20','-2','4','-4','2');

  default = str2mat('default',default1,default2);
  
  harmonic1 = str2mat(...
      'x','v',...
      ' v',' -(k*x + d*v)/m',...
      'k','d','m',' ');
  harmonic2 = str2mat(...
      '3','0','1',' ',...
      '2','20','-5','5','-5','5');
  
  harmonic = str2mat('harmonic',harmonic1, harmonic2);

  pend1 = str2mat(...
      'theta','ang_vel',...
      ' ang_vel',' - sin(theta) -D*ang_vel',...
      'D',' ',' ','  ');
  pend2 = str2mat(...
      '0',' ',' ',' ',...
      '2','20','-10','10','-4','4');
  
  pendulum = str2mat('pendulum',pend1,pend2);
	
  predprey1 = str2mat(...
      'prey','predator',...
      ' A*prey - B*prey*predator',' -C*predator + D*prey*predator',...
      'A','B','C','D');
  predprey2 = str2mat(...
      '0.4','0.01','0.3','0.005',...
      '2','20','0','120','0','80');
  predprey = str2mat('prepre',predprey1,predprey2);
	
  linear1 = str2mat(...
      'x','y',...
      ' A*x + B*y',' C*x + D*y',...
      'A','C','B','D');
  linear2 =str2mat(...
      '2','-2','2','-3',...
      '2','20','-5','5','-5','5');
  
  linear = str2mat('linear',linear1,linear2);
	
  compspec1 = str2mat('x','y',...
      ' (1 - x - y)*x',' (A - B*x - C*y)*y',...
      'A','B','C',' ');
  compspec2 = str2mat('4','2','7',' ',...
      '2','20','0','1','0','1');
  
  compspec = str2mat('compspec',compspec1,compspec2);
	
  vanderpol1 = str2mat('x','y',...
      ' M*x - y -x^3',' x',...
      'M',' ',' ',' ');
  vanderpol2 = str2mat('2','0','0','0',...
      '2','20','-5','5','-5','5');
  
  vanderpol = str2mat('vanderpol', vanderpol1,vanderpol2);
	
  duffing1 = str2mat('x','y',...
      ' y',' -(k*x + c*y + l*x^3)/m',...
      'k','m','c','l');
  duffing2 = str2mat('-1','1','0','1',...
      '2','20','-2','2','-2','2');
  duffing = str2mat('duff',duffing1,duffing2);
		
  systems = str2mat(default,harmonic,pendulum,...
      predprey,linear,compspec,vanderpol,duffing);

  

  
  NumbFPts = 20;	% The number of field points plotted in each direction.
  Arrflag = 2;		% The default is a vector field.
  magn = 1.25;		% The expansion factor in each direction 
                        % for the compute window.
  steps = 20;
  tol = 1e-6;
  itlim = 1300;
  settings = [Arrflag,NumbFPts,magn,steps,itlim,tol];
  WINvect = [-4 4 -4 4];
  
  Xname = ' x';
  Yname = ' y';
  xderivstr = ' 2*x - y + 3*(x^2-y^2) + 2*x*y';
  yderivstr = ' x - 3*y - 3*(x^2-y^2) + 3*x*y';
  WINvect=[-2 4 -4 2];
  pname1='';
  pname2='';
  pname3='';
  pname4='';
  parav1='';
  parav2='';
  parav3='';
  parav4='';
  
  texth =22;   		% Height of text boxes.
  varw = 80;		% Length of variable boxes.
  equalw =30;		% Length of equals.
  eqlength = 350;	% Length of right hand sides of equations.
  winstrlen = 180;	% Length of string boxes in display frame.
  left =25;		% Left margin of the frames.
  frsep = 10;		% Separation between frames.
  separation = texth+3;	% Separation between boxes.
  
  dfigwidth =2*left + varw+equalw+eqlength+10;	% Width of the figure.
  dfigurebot = 30;				% Bottom of the figure.
  
  qwind = [dfigwidth/4-30,frsep,60,separation];	% Quit button
  rwind = [dfigwidth/2-30,frsep,60,separation];	% Revert "
  pwind = [3*dfigwidth/4-30,frsep,60,separation];	% Proceed "
  
  disfrbot = 2*frsep + separation;	% Display frame.
  disfrw = winstrlen + 50 +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;
  ffrleft = left + disfrw + frsep;
  ffrw = dfigwidth -left - ffrleft;
  ffrht = disfrht; 
  ffrwind = [ffrleft, ffrbot, ffrw, ffrht];
  
  dfigureheight = defrbot + defrht +frsep;	% Height of the figure.
  
  ppset = figure('pos',[30 dfigurebot dfigwidth dfigureheight],...
      'resize','off','name','PPLANE Setup','numb','off','visible','off');
  figure(ppset)
  
  
  frame(1) = uicontrol('style','frame','pos',defrwind,'visible','off');
  
  xname=[
    'sh = get(gcf,''user'');'...
	'sm = get(sh(1),''user'');'...
	'eq = sh(4:10);'...
	'wind =sh(11:19);'...
	'Xname=get(eq(2),''string'');'...
	'set(wind(2),''string'',[''The minimum value of '',Xname,'' = '']);'...
	'set(wind(4),''string'',[''The maximum value of '',Xname,'' = '']);'...
	'sm = str2mat(Xname,sm(2:18,:));'...
	'set(sh(1),''user'',sm);'...
	'set(sh(3),''user'',0);'];
  
  yname=[
    'sh = get(gcf,''user'');'...
	'sm = get(sh(1),''user'');'...
	'eq = sh(4:10);'...
	'wind =sh(11:19);'...
	'Yname=get(eq(5),''string'');'...
	'set(wind(6),''string'',[''The minimum value of '',Yname,'' = '']);'...
	'set(wind(8),''string'',[''The maximum value of '',Yname,'' = '']);'...
	'sm = str2mat(sm(1,:),Yname,sm(3:18,:));'...
	'set(sh(1),''user'',sm);'...
	'set(sh(3),''user'',0);'];
  
  xder =[
    'sh = get(gcf,''user'');'...
	'sm = get(sh(1),''user'');'...
	'eq=sh(4:10);'...
	'xderivstr = get(eq(4),''string'');'...
	'sm = str2mat(sm(1:2,:),xderivstr,sm(4:18,:));'...
	'set(sh(1),''user'',sm);'...
	'set(sh(3),''user'',0);'];
  
  yder =[
    'sh = get(gcf,''user'');'...
	'sm = get(sh(1),''user'');'...
	'eq=sh(4:10);'...
	'yderivstr = get(eq(7),''string'');'...
	'sm = str2mat(sm(1:3,:),yderivstr,sm(5:18,:));'...
	'set(sh(1),''user'',sm);'...
	'set(sh(3),''user'',0);'];
  
  equationbot = defrbot + 5;
  eqlabelbot = equationbot + 2*separation;
  xbot = equationbot + separation;		% Bottom of x equation.
  ybot = equationbot;				% Bottom of y equation.
  
  lablen =300;
  
  eqlableft = (dfigwidth-lablen)/2;
  eqleft = left + 5;
  
  eq(1)=uicontrol('style','text',...
      'pos',[eqlableft eqlabelbot lablen texth],...
      'horizon','center',...
      'string','The differential equations.','visible','off');
  
  eq(2)=uicontrol('pos',[eqleft, xbot, varw, texth],'style','edit',...
      'horizon','right','string',Xname,'call',xname,'visible','off');
  
  eq(3) = uicontrol('style','text','pos',[eqleft+varw xbot equalw texth],...
      'horizon','center','string',''' = ','visible','off');
  
  eq(4)=uicontrol('pos',[eqleft+varw+equalw xbot eqlength texth],...
      'string',xderivstr,...
      'horizon','left','style','edit',...
      'call',xder,'visible','off');
  
  eq(5)=uicontrol('pos',[eqleft ybot varw texth],'style','edit',...
      'horizon','right','string',Yname,'call',yname,'visible','off');
  
  eq(6) = uicontrol('style','text',...
      'pos',[eqleft+varw ybot equalw texth],...
      'horizon','center','string',''' = ','visible','off');
  
  eq(7)=uicontrol('pos',[eqleft+varw + equalw ybot eqlength texth],...
      'string',yderivstr,...
      'horizon','left','style','edit',...
      'call',yder,'visible','off');
  
  
  frame(2) = uicontrol('style','frame','pos',disfrwind,'visible','off');
  
  w1 = [
    'sh = get(gcf,''user'');'...
	'wind =sh(11:19);'...
	'nnn = str2num(get(wind(3),''string''));',...
	'if isempty(nnn),',...
	'  set(wind(3),''string'',''?'');',...
	'  nnn = NaN;',...
	'end,',...
	'sh(36)=nnn;'...
	'set(gcf,''user'',sh);',...
	'sm = get(sh(1),''user'');',...
	'sm = str2mat(sm(1:14,:),get(wind(3),''string''),sm(16:18,:));',...
	'set(sh(1),''user'',sm);'];	
  
  w2 = [
    'sh = get(gcf,''user'');'...
	'wind =sh(11:19);'...
	'nnn = str2num(get(wind(5),''string''));',...
	'if isempty(nnn),',...
	'  set(wind(5),''string'',''?'');',...
	'  nnn = NaN;',...
	'end,',...
	'sh(37)=nnn;'...
	'set(gcf,''user'',sh);',...
	'sm = get(sh(1),''user'');',...
	'sm = str2mat(sm(1:15,:),get(wind(5),''string''),sm(17:18,:));',...
	'set(sh(1),''user'',sm);'];	
  
  w3 = [
    'sh = get(gcf,''user'');'...
	'wind =sh(11:19);'...
	'nnn = str2num(get(wind(7),''string''));',...
	'if isempty(nnn),',...
	'  set(wind(7),''string'',''?'');',...
	'  nnn = NaN;',...
	'end,',...
	'sh(38)=nnn;',...
	'set(gcf,''user'',sh);',...
	'sm = get(sh(1),''user'');',...
	'sm = str2mat(sm(1:15,:),get(wind(5),''string''),sm(17:18,:));',...
	'set(sh(1),''user'',sm);'];	

  w4 = [
    'sh = get(gcf,''user'');'...
	'wind =sh(11:19);'...
	'nnn = str2num(get(wind(9),''string''));',...
	'if isempty(nnn),',...
	'  set(wind(9),''string'',''?'');',...
	'  nnn = NaN;',...
	'end,',...
	'sh(39)=nnn;'...
	'set(gcf,''user'',sh);',...
	'sm = get(sh(1),''user'');',...
	'sm = str2mat(sm(1:16,:),get(wind(7),''string''),sm(18,:));',...
	'set(sh(1),''user'',sm);'];	
	
  winbot1 = disfrbot + disfrht - 5 - separation;
  winbot2 = winbot1 - separation;
  winbot3 = winbot2 - separation;
  winbot4 = winbot3 - separation;
  winbot5 = winbot4 - separation;
  
  wind(1)=uicontrol('style','text',...
      'pos',[eqleft winbot1 disfrw-10 texth],...
      'horizon','center',...
      'string','The display window.','visible','off');
  
  wind(2) = uicontrol('style','text',...
      'pos',[eqleft winbot2 winstrlen texth],...
      'horizon','right',...
      'string',['The minimum value of ',Xname,' = '],'visible','off');
  
  wind(3) = uicontrol('style','edit',...
      'pos',[eqleft+winstrlen winbot2 50 texth],...
      'string',num2str(WINvect(1)),...
      'call',w1,'visible','off');
  
  wind(4) = uicontrol('style','text',...
      'pos',[eqleft winbot3 winstrlen texth],...
      'horizon','right',...
      'string',['The maximum value of ',Xname,' = '],'visible','off');
  
  wind(5) = uicontrol('style','edit',...
      'pos',[eqleft+winstrlen winbot3 50 texth],...
      'string',num2str(WINvect(2)),...
      'call',w2,'visible','off');
  
  wind(6) = uicontrol('style','text',...
      'pos',[eqleft winbot4 winstrlen texth],...
      'horizon','right',...
      'string',['The minimum value of ',Yname,' = '],'visible','off');

  wind(7) = uicontrol('style','edit',...
      'pos',[eqleft+winstrlen winbot4 50 texth],...
      'string',num2str(WINvect(3)),...
      'call',w3,'visible','off');
	
  wind(8) = uicontrol('style','text',...
      'pos',[eqleft winbot5 winstrlen texth],...
      'horizon','right',...
      'string',['The maximum value of ',Yname,' = '],'visible','off');
		
  wind(9) = uicontrol('style','edit',...
      'pos',[eqleft+winstrlen winbot5 50 texth],...
      'string',num2str(WINvect(4)),...
      'call',w4,'visible','off');
  
  
  frame(3)=uicontrol('style','frame','pos',pfrwind,'visible','off');
  
  p1 = [
    'sh = get(gcf,''user'');'...
	'sm = get(sh(1),''user'');'...
	'eq=sh(4:10);'...
	'wind =sh(11:19);'...
	'frame=sh(1:3);'...
	'para=sh(20:32);'...
	'pname1=get(para(2),''string'');'...
	'parav1=get(para(4),''string'');'...
	'sm = str2mat(sm(1:4,:),pname1,sm(6:8,:),parav1,sm(10:18,:));'...
	'set(sh(1),''user'',sm);'...
	'set(sh(3),''user'',0);'];
  
  p2 = [
    'sh = get(gcf,''user'');'...
	'sm = get(sh(1),''user'');'...
	'eq=sh(4:10);'...
	'wind =sh(11:19);'...
	'frame=sh(1:3);'...
	'para=sh(20:32);'...
	'pname2=get(para(5),''string'');'...
	'parav2=get(para(7),''string'');'...
	'sm = str2mat(sm(1:5,:),pname2,sm(7:9,:),parav2,sm(11:18,:));'...
	'set(sh(1),''user'',sm);'...
	'set(sh(3),''user'',0);'];

  p3 = [
    'sh = get(gcf,''user'');'...
	'sm = get(sh(1),''user'');'...
	'eq=sh(4:10);'...
	'wind =sh(11:19);'...
	'frame=sh(1:3);'...
	'para=sh(20:32);'...
	'pname3=get(para(8),''string'');'...
	'parav3=get(para(10),''string'');'...
	'sm = str2mat(sm(1:6,:),pname3,sm(8:10,:),parav3,sm(12:18,:));'...
	'set(sh(1),''user'',sm);'...
	'set(sh(3),''user'',0);'];

  p4 = [
    'sh = get(gcf,''user'');'...
	'sm = get(sh(1),''user'');'...
	'eq=sh(4:10);'...
	'wind =sh(11:19);'...
	'frame=sh(1:3);'...
	'para=sh(20:32);'...
	'pname4=get(para(11),''string'');'...
	'parav4=get(para(13),''string'');'...
	'sm = str2mat(sm(1:7,:),pname4,sm(9:11,:),parav4,sm(13:18,:));'...
	'set(sh(1),''user'',sm);'...
	'set(sh(3),''user'',0);'];

  pvarw = 75;
  peqw = 15;
  pbot2 = pfrbot + 5;
  pbot1 = pbot2 + separation;
  pleft1 = eqleft + 80 + 10;
  peqleft1 = pleft1 + pvarw;
  pvleft1 = peqleft1 + peqw;
  pleft2 = pvleft1 + 50 +10;
  peqleft2 = pleft2 + pvarw;
  pvleft2 = peqleft2 + peqw;	
  
  para(1)=uicontrol('style','text',...
      'pos',[eqleft pbot2+separation/2 80 texth],...
      'horizon','center',...
      'string','Parameters:','visible','off');
  
  para(2) = uicontrol('style','edit',...
      'pos',[pleft1 pbot1 pvarw texth],...
      'horizon','right','string',pname1,...
      'call',p1,'visible','off');
  
  para(3) = uicontrol('style','text',...
      'pos',[peqleft1 pbot1 peqw texth],...
      'horizon','center','string','=','visible','off');
  
  para(4) = uicontrol('style','edit',...
      'pos',[pvleft1 pbot1 50 texth],...
      'string', '','call',p1,'visible','off');
  
  para(5) = uicontrol('style','edit',...
      'pos',[pleft1 pbot2 pvarw texth],...
      'horizon','right','string',pname2,...
      'call',p2,'visible','off');
  
  para(6) = uicontrol('style','text',...
      'pos',[peqleft1 pbot2 peqw texth],...
      'horizon','center','string','=','visible','off');
  
  para(7) = uicontrol('style','edit',...
      'pos',[pvleft1 pbot2 50 texth],...
      'string','',...
      'call',p2,'visible','off');
  
  para(8) = uicontrol('style','edit',...
      'pos',[pleft2 pbot1 pvarw texth],...
      'horizon','right','string',pname3,...
      'call',p3,'visible','off');
  
  para(9) = uicontrol('style','text',...
      'pos',[peqleft2 pbot1 peqw texth],...
      'horizon','center','string','=','visible','off');
  
  para(10) = uicontrol('style','edit',...
      'pos',[pvleft2 pbot1 50 texth],...
      'string','','call',p3,'visible','off');
  
  para(11) = uicontrol('style','edit',...
      'pos',[pleft2 pbot2 pvarw texth],...
      'horizon','right','string',pname4,...
      'call',p4,'visible','off');
  
  para(12) = uicontrol('style','text',...
      'pos',[peqleft2 pbot2 peqw texth],...
      'horizon','center','string','=','visible','off');
  
  para(13) = uicontrol('style','edit',...
      'pos',[pvleft2 pbot2 50 texth],...
      'string','','call',p4,'visible','off');
  
  butt(1) = uicontrol('style','push',...
      'pos',qwind,...
      'string','Quit','call','pplane(''quit'')','visible','off');
  
  butt(2) = uicontrol('style','push',...
      'pos',rwind,...
      'string','Revert','call','pplane(''revert'')','visible','off');
	    
  butt(3) = uicontrol('style','push',...
      'pos',pwind,...
      'string','Proceed','call','pplane(''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 + 10;
  radw = 80;
  
  textleft = radleft + radw + 15;
  textw = ffrleft + ffrw -textleft -5;
  
  typewind = [ffrleft+5, ffrbot, radw+10, ffrht-separation-7];
  
  textwind = [textleft-5, ffrbot,textw+5, ffrht-separation-7];
  
  typeframe = uicontrol('style','frame','pos',typewind,'visible','off');
  
  textframe = uicontrol('style','frame','pos',textwind,'visible','off');
  
  
  
  rval1 = ~(Arrflag -1);
  rval2 = (~(Arrflag -2))*2;
  rval3 = (~(Arrflag -3))*3;
  
  rtitle = uicontrol('style','text',...
      'pos',[radleft winbot2 radw texth],...
      'string','Type.','horizon','center','visible','off');
  
  rad(1) = uicontrol('style','radio',...
      'pos',[radleft winbot3 radw texth],...
      'string','Lines','value',rval1,'visible','off');
  
  rad(2) = uicontrol('style','radio',...
      'pos',[radleft winbot4 radw texth],...
      'string','Arrows','value',rval2,...
      'max',2,'visible','off');

  rad(3) = uicontrol('style','radio',...
      'pos',[radleft winbot5 radw texth],...
      'string','None','value',3*rval3,'max',3,'visible','off');
				
  for i=1:3
    set(rad(i),'user',rad(:,[1:(i-1),(i+1):3]));
  end
  
  callrad = [
    'me = get(gcf,''currentobject'');',...
	'kk = get(me,''max'');',...
	'set(get(me,''user''),''value'',0),',...
	'set(me,''value'',kk),',...
	'sh = get(gcf,''user'');',...
	'sh(44) = kk;',...
	'set(gcf,''user'',sh);',...
	'sm = get(sh(1),''user'');',...
	'sm = str2mat(sm(1:12,:),num2str(kk),sm(14:18,:));',...
	'set(sh(1),''user'',sm);'];
  
  set(rad,'call',callrad);
  
  text1 = uicontrol('style','text',...
      'pos',[textleft winbot2 textw texth],...
      'string','Number of','horizon','left','visible','off');
  
  text2 = uicontrol('style','text',...
      'pos',[textleft winbot3 textw texth],...
      'string','field points per','horizon','left','visible','off');

  text3 = uicontrol('style','text',...
      'pos',[textleft winbot4 textw texth],...
      'string','row or column.','horizon','left','visible','off');

  callnfpts = [
    'me = get(gcf,''currentobject'');',...
	'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,',...
	'sh = get(gcf,''user'');'...
	'sh(45) = kk;',...
	'set(gcf,''user'',sh);'];
  
  nfpts = uicontrol('style','edit',...
      'pos',[textleft+(textw -50)/2,winbot5 50,texth],...
      'string',num2str(NumbFPts),'call',callnfpts,'visible','off');

  mesyst = uimenu(ppset,'label','Systems','vis','off');
  meload = uimenu(mesyst,'label','Load system',...
      'call','pplane(''loadsyst'');');
  mesave = uimenu(mesyst,'label','Save system',...
      'call','pplane(''savesyst'');');
  meadd = uimenu(mesyst,'label','Add current system to the gallery',...
      'call','pplane(''addgall'');');
  
  sysmenu = uimenu('label','Gallery',...
      'visible','off','tag','gallery');
  
  medef = uimenu(sysmenu,'label','default system',...
      'call','pplane(''system'',''default'')','visible','off');
  melin = uimenu(sysmenu,'label','linear system',...
      'call','pplane(''system'',''linear'')','visible','off');
  meharm = uimenu(sysmenu,'label','vibrating spring  ',...
      'call','pplane(''system'',''harmonic'')','visible','off');
  mepend = uimenu(sysmenu,'label','pendulum',...
      'call','pplane(''system'',''pendulum'')','visible','off');
  meprepre = uimenu(sysmenu,'label','predator prey',...
      'call','pplane(''system'',''prepre'')','visible','off');
  mecompspec = uimenu(sysmenu,'label','competing species',...
      'call','pplane(''system'',''compspec'')','visible','off');
  mevdp = uimenu(sysmenu,'label','van der Pol''s equation',...
      'call','pplane(''system'',''vanderpol'')','visible','off');
  meduff = uimenu(sysmenu,'label','Duffing''s equation',...
      'call','pplane(''system'',''duff'')','visible','off');
  
  % Record the handles in the User Data of the Set Up figure.
  
  sethandles = [frame,eq,wind,para,butt,WINvect,WINvect,settings,rad,nfpts];
  set(ppset,'user',sethandles);
  hhhh = get(ppset,'child');
  
  % uimenus do not have the 'units' property
  
  kkk = find (hhhh == sysmenu);       
  hhhh = findobj(ppset,'type','uicontrol');
  set(hhhh,'units','normal')
  
  global PPSETPOS
  
  if all(size(PPSETPOS) == [1,4])
    set(gcf,'pos',PPSETPOS);
  end
  
  set(ppset,'visible','on','resize','on');
  set(get(ppset,'children'),'visible','on');
  set(get(sysmenu,'children'),'visible','on');
  
  % Record the important information in the User Data of
  % varous controls.  This will be dynamic storage.
  
  af = num2str(Arrflag);
  npts = num2str(NumbFPts);
  aaa = num2str(WINvect(1));
  bbb = num2str(WINvect(2));
  ccc = num2str(WINvect(3));
  ddd = num2str(WINvect(4));
  
  sm = str2mat(Xname,Yname,xderivstr,yderivstr,pname1,pname2,...
      pname3,pname4,str2mat(parav1,parav2,parav3,parav4,af,...
      npts,aaa,bbb,ccc,ddd));
  set(sethandles(1),'user',sm);
  set(sethandles(2),'user',sm);
  set(sethandles(3),'user',0);
  set(sethandles(4),'user',systems);
  
  eval(get(eq(4),'call'))
  eval(get(eq(7),'call'))

elseif strcmp(action,'savesyst')
  
  sh = get(gcf,'user');
  systems = get(sh(4),'user');
  [fname,pname] = uiputfile('*.sys','Save as:');
  if (strcmp(fname,'') | fname == 0)
    return
  end    
  ll = length(fname);
  if (ll<=4)
    fn = fname;    
  elseif (ll>4 & strcmp(fname(ll-3:ll),'.sys'))
    fn = fname(1:(ll-4));
  else
    fn =fname;
  end
  newsyst = get(sh(1),'user');
  newsyst = str2mat(fn,newsyst);
  ll = size(newsyst,1);
  fid = fopen([pname fname],'w');
  for kk = 1:ll
    fprintf(fid,[newsyst(kk,:),'\n']);
  end
  fclose(fid);

elseif strcmp(action,'loadsyst')
  
  sh = get(gcf,'user');
  systems = get(sh(4),'user');
  [fname,pname] = uigetfile('*.sys','Select a system to load.');

  if (strcmp(fname,'') | fname == 0)
    return
  end    
  if isstr(fname)
    fid = fopen([pname fname],'r');
    newsyst = fgetl(fid);
    while ~feof(fid)
      newsyst = str2mat(newsyst, fgetl(fid));
    end
    kk = size(newsyst,1);
    if kk >= 19
      newsyst = newsyst([1:19],:);

      systems = str2mat(systems,newsyst);
      set(sh(4),'user',systems);
      pplane('system',deblank(newsyst(1,:)));
      pplane('addgall',newsyst);
    else
      warndlg(['The file ',fname, 'does define a proper system.']);
    end
  end

elseif strcmp(action,'addgall')

  if nargin == 2
    sh = get(gcf,'user');
    systems = get(sh(4),'user');
    gallery = findobj('tag','gallery');
    syst = input1;
    sname = deblank(syst(1,:));
    if length(get(gallery,'child'))<9
      sepstr = 'on';
    else
      sepstr = 'off';
    end
    newmenu = uimenu(gallery,'label',sname,...
      'call',['pplane(''system'',''',sname,''');'],...
      'separator',sepstr);

  else  
    edge = 5; fudge = 3; sep = 5;
    sw = 300;
    ht = 20;
    sl = edge + fudge;
    bb = sep + edge;
    nb = bb + ht + sep;
    pb = nb + ht + sep;
    frw = sw + 2*fudge;
    figw = frw + 2*edge;
    frh = 4*sep + 3*ht;
    figh = frh + 2*edge;
    figpos = [200,300,figw,figh];
    frpos = [edge,edge,frw,frh];
    buttw = sw/2;
    canpos = [sl,bb,buttw,ht];
    addpos = [sl+buttw,bb,buttw,ht];
    npos = [sl,nb,sw,ht];
    ppos = [sl,pb,sw,ht];
    
    dlg = figure('pos',figpos,'numb','off',...
	'name','PPLANE Add text',...
	'tag','PPLANE','resize','off',...
	'vis','off');
    
    dfr = uicontrol('style','frame','pos',frpos);
    
    prompt = uicontrol('style','text','pos',ppos,...
	'string','Give the system a name.',...
	'horiz','left');
    
    
    name = uicontrol('style','edit','pos',npos,...
	'horiz','left','tag','systemname');
    
    cancel = uicontrol('style','push','pos',canpos,...
	'string','Cancel','call','delete(gcf)');
    
    addcall = [
      'me = findobj(''tag'',''systemname'');',...
	  'mee = gcf;',...
	  'ppset = findobj(''name'',''PPLANE Setup'');',...
	  'sname = get(me,''string'');',...
	  'sh = get(ppset,''user'');',...
	  'systems = get(sh(4),''user'');',...
	  'gallery = findobj(''tag'',''gallery'');',...
	  'syst = get(sh(1),''user'');',...
	  'syst = str2mat(sname,syst);',...
	  'systems = str2mat(systems,syst);',...
	  'set(sh(4),''user'',systems);',...
	  'if length(get(gallery,''child''))<9,',...
	  '  sepstr = ''on'';',...
	  'else,',...
	  '  sepstr = ''off'';',...
	  'end,',...
	  'newmenu = uimenu(gallery,''label'',sname,',...
	  '  ''call'',[''pplane(''''system'''','''''',sname,'''''');''],',...
	  '  ''separator'',sepstr);',...
	  'delete(mee)'];
    
    add = uicontrol('style','push','pos',addpos,...
	'string','Add','call',addcall);
    
    h = [dlg,dfr,prompt,name,cancel,add];
    set(h,'units','normal','vis','on');
    set(dlg,'resize','on');
  end

elseif strcmp(action,'system')
  
  sh = get(gcf,'user');
  set(sh(3),'user',0);
  systems = get(sh(4),'user');
  ll = size(systems,1);
  if rem(ll,19)
    error('There is something wrong with the list of systems.')
  end
  k=1;
  while ~strcmp(input1,deblank(systems(19*(k-1)+1,:)))
    k = k+1;
  end
  syst = systems(19*(k-1)+2:19*k,:);


  frame=sh(1:3);
  eq=sh(4:10);
  wind =sh(11:19);
  para=sh(20:32);
  
  Xname = deblank(syst(1,:));
  Yname = deblank(syst(2,:));
  xderivstr = deblank(syst(3,:));
  yderivstr = deblank(syst(4,:));
  pname1 = deblank(syst(5,:));
  pname2 = deblank(syst(6,:));
  pname3 = deblank(syst(7,:));
  pname4 = deblank(syst(8,:));
  parav1 = deblank(syst(9,:));
  parav2 = deblank(syst(10,:));
  parav3 = deblank(syst(11,:));
  parav4 = deblank(syst(12,:));
  Arrflag = str2num(syst(13,:));
  NumbFPts = str2num(syst(14,:));
  a = str2num(syst(15,:));
  b = str2num(syst(16,:));	
  c = str2num(syst(17,:));
  d = str2num(syst(18,:));
  
  WINvect = [a,b,c,d];
  sm = syst;
  
  set(sh(1),'user',sm);
  
  set(eq(2),'string',Xname);
  set(wind(2),'string',['The minimum value of ',Xname,' = ']);
  set(wind(4),'string',['The maximum value of ',Xname,' = ']);
  set(eq(5),'string',Yname);
  set(wind(6),'string',['The minimum value of ',Yname,' = ']);
  set(wind(8),'string',['The maximum value of ',Yname,' = ']);
  set(eq(4),'string',xderivstr);
  set(eq(7),'string',yderivstr);
  set(para(2),'string',pname1);
  set(para(5),'string',pname2);
  set(para(8),'string',pname3);
  set(para(11),'string',pname4);
  set(para(4),'string',parav1);
  set(para(7),'string',parav2);
  set(para(10),'string',parav3);
  set(para(13),'string',parav4);
  set(wind(3),'string',num2str(a));
  set(wind(5),'string',num2str(b));
  set(wind(7),'string',num2str(c));
  set(wind(9),'string',num2str(d));
  sh(36:39) = WINvect;
  sh(44) = Arrflag;
  sh(45) = NumbFPts;
  set(gcf,'user',sh);
  
  rval = zeros(1,3);
  rval(1) = ~(Arrflag -1);
  rval(2) = (~(Arrflag -2))*2;
  rval(3) = (~(Arrflag -3))*3;
  
  rad = sh(50:52);
  
  for i=1:3
    set(rad(i),'value',rval(i));
  end
  
  nfpts = sh(53);
  set(nfpts,'string',deblank(syst(14,:)));
	

elseif strcmp(action,'proceed')

  % Proceed connects Setup with the Display window.

  ppset = gcf;
  sh = get(ppset,'user');
  flag = get(sh(3),'user');
  
  %  flag = 0, if this is the first time through for this equation, 
  %  but flag = 1 if only the window dimensions or the field data 
  %  have been changed.

  WINvect = sh(36:39);
  if any(isnan(WINvect))
    errmsg = ['One of the entries defining the display window ',...
	  'is not a number.'];
    fprintf('\a')
    errordlg(errmsg,'PPLANE error','on');
    return
  end
  sm = get(sh(1),'user');
  xstr = deblank(sm(1,:));
  if isempty(xstr)
    errmsg = 'The first dependent variable needs a name.';
    fprintf('\a')
    errordlg(errmsg,'PPLANE error','on');
    return
  end
  ystr = deblank(sm(2,:));
  if isempty(ystr)
    errmsg = 'The second dependent variable needs a name.';
    fprintf('\a')
    errordlg(errmsg,'PPLANE error','on');
    return
  end
  if WINvect(2)<= WINvect(1)
    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)
    errmsg = ['The minimum value of ', ystr,...
	  ' must be smaller than the maximum value.'];
    fprintf('\a')
    errordlg(errmsg,'PPLANE error','on');
    return
  end
  XxXxXx = (WINvect(1)+WINvect(2))/2;
  YyYyYy = (WINvect(3)+WINvect(4))/2;
  err = 0;
  eval([xstr,'=XxXxXx;'],'err = 1;');
  if err
    errmsg = ['"',xstr, '" is not a valid variable name in MATLAB.'];
    fprintf('\a')
    errordlg(errmsg,'PPLANE error','on');
    return
  end
  err = 0;
  eval([ystr,'=YyYyYy;'],'err = 1;');
  if err
    errmsg = ['"',ystr, '" is not a valid variable name in MATLAB.'];
    fprintf('\a')
    errordlg(errmsg,'PPLANE error','on');
    return
  end
  Arrflag = sh(44);
  NumbFPts = sh(45);
  if isnan(NumbFPts)
    errmsg = 'The entry for the number of field points is not a number.';
    fprintf('\a')
    errordlg(errmsg,'PPLANE error','on');
    return
  end
    
  
  if (flag == 1)
    ppdisp = findobj('name','PPLANE Display');
    if isempty(ppdisp)   %  this is a real kluge.
      set(sh(3),'user',1);
      set(ppset,'user',sh);	
      pplane('display');
      break;
    end      
    dh = get(ppdisp,'user');
    
    if ( (NumbFPts ~= dh(24)) | (any(WINvect ~= dh(19:22)) ) )
      sh(40:43) = WINvect;			
      dh(19:22) = WINvect;
      dh(23:24) = [Arrflag, NumbFPts];
      set(ppdisp,'user',dh);
      pplane('dirfield');
      
    elseif (Arrflag ~= dh(23))
      arrh = get(dh(6),'user');
      if (Arrflag == 1)
	set(arrh(1),'visible','on');
	set(arrh(2),'visible','off');
      elseif(Arrflag == 2)
	set(arrh(1),'visible','off');
	set(arrh(2),'visible','on');
      else
	set(arrh,'visible','off');
      end
      dh(23) = Arrflag;
      set(ppdisp,'user',dh);
      pplane('refresh')
    end
  else
    set(sh(3),'user',1);
    set(ppset,'user',sh);	
    pplane('display');
  end	
		
elseif strcmp(action,'revert')

	sh = get(gcf,'user');
	frame=sh(1:3);
	eq=sh(4:10);
	wind =sh(11:19);
	para=sh(20:32);
	sm = get(sh(2),'user');
	Xname = deblank(sm(1,:));
	Yname = deblank(sm(2,:));
	xderivstr = deblank(sm(3,:));
	yderivstr = deblank(sm(4,:));
	pname1 = deblank(sm(5,:));
	pname2 = deblank(sm(6,:));
	pname3 = deblank(sm(7,:));
	pname4 = deblank(sm(8,:));
	parav1 = deblank(sm(9,:));
	parav2 = deblank(sm(10,:));
	parav3 = deblank(sm(11,:));
	parav4 = deblank(sm(12,:));
	Arrflag = str2num(sm(13,:));
	NumbFPts = str2num(sm(14,:));
	WINvect = sh(40:43);

	set(sh(1),'user',sm);
	
	set(eq(2),'string',Xname);
	set(wind(2),'string',['The minimum value of ',Xname,' = ']);
	set(wind(4),'string',['The maximum value of ',Xname,' = ']);
	set(eq(5),'string',Yname);
	set(wind(6),'string',['The minimum value of ',Yname,' = ']);
	set(wind(8),'string',['The maximum value of ',Yname,' = ']);
	set(eq(4),'string',xderivstr);
	set(eq(7),'string',yderivstr);
	set(para(2),'string',pname1);
	set(para(5),'string',pname2);
	set(para(8),'string',pname3);
	set(para(11),'string',pname4);
	set(para(4),'string',parav1);
	set(para(7),'string',parav2);
	set(para(10),'string',parav3);
	set(para(13),'string',parav4);
	set(wind(3),'string',num2str(WINvect(1)));
	set(wind(5),'string',num2str(WINvect(2)));
	set(wind(7),'string',num2str(WINvect(3)));
	set(wind(9),'string',num2str(WINvect(4)));
	sh(36:39) = WINvect;
	sh(44) = Arrflag;
	sh(45) = NumbFPts;
	set(gcf,'user',sh);

	rval = zeros(1,3);
	rval(1) = ~(Arrflag -1);
	rval(2) = (~(Arrflag -2))*2;
	rval(3) = (~(Arrflag -3))*3;

	rad = sh(50:52);
				
	for i=1:3
		set(rad(i),'value',rval(i));
	end
	
	nfpts = sh(53);
	set(nfpts,'string',deblank(sm(14,:)));
	

elseif strcmp(action,'display')

  % Initialize display window.

  % Display takes the information from the Setup window and initializes
  % the Display window if it already exists.  If the Display window does
  % not exist, it builds one.

  ppset = gcf;

  % Get the information from PPLANE Setup

  sh = get(ppset,'user');
  settings = sh(44:49);
  WINvect = sh(36:39);
	
  sm = get(sh(1),'user');
  set(sh(2),'user',sm);	% Resets long term storage.
	
  Xname = deblank(sm(1,:));
  Yname = deblank(sm(2,:));
  xderivstr = deblank(sm(3,:));
  yderivstr = deblank(sm(4,:));
  pname1 = deblank(sm(5,:));
  pname2 = deblank(sm(6,:));
  pname3 = deblank(sm(7,:));
  pname4 = deblank(sm(8,:));
  parav1 = deblank(sm(9,:));
  parav2 = deblank(sm(10,:));
  parav3 = deblank(sm(11,:));
  parav4 = deblank(sm(12,:));
  sh(40:43) = WINvect;
  set(ppset,'user',sh);
  
  % 	Convert the parameters to their current values.
  
  %	First remove the blanks.
  xderivstr(find(abs(xderivstr)==32))=[];
  yderivstr(find(abs(yderivstr)==32))=[];
  
  %       Next remove the periods inserted by users attempting to make 
  %       the function array smart.
  
  l=length(xderivstr);
  for ( k = fliplr(findstr('.',xderivstr)))
    if k == l
      xderivstr = xderivstr(1:l-1);
    elseif (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 k == l
      yderivstr = yderivstr(1:l-1)
    elseif (find('*/^' == yderivstr(k+1)))
      yderivstr = [yderivstr(1:k-1), yderivstr(k+1:l)];
    end
    l=l-1;
  end
  
  %	Replace the parameters with their values.
  perr = [];
  if (pname1 ~= [])
    if isempty(str2num(parav1))
      perr = 23;
    else
      xderivstr = pplane('paraeval',pname1,parav1,xderivstr);
      yderivstr = pplane('paraeval',pname1,parav1,yderivstr);
    end
  end
  if (pname2 ~= [])
    if isempty(str2num(parav2))
      perr = [perr;26];
    else
      xderivstr = pplane('paraeval',pname2,parav2,xderivstr);
      yderivstr = pplane('paraeval',pname2,parav2,yderivstr);
    end
  end
  if (pname3 ~= [])
    if isempty(str2num(parav3))
      perr = [perr;29];
    else
      xderivstr = pplane('paraeval',pname3,parav3,xderivstr);
      yderivstr = pplane('paraeval',pname3,parav3,yderivstr);
    end
  end
  if( pname4 ~= [])
    if isempty(str2num(parav4))
      perr = [perr;32];
    else
      xderivstr = pplane('paraeval',pname4,parav4,xderivstr);
      yderivstr = pplane('paraeval',pname4,parav4,yderivstr);
    end
  end

  %  Do the derivative strings make sense?

  XxXxXx = WINvect(1)+rand*(WINvect(2) - WINvect(1));
  YyYyYy = WINvect(3)+rand*(WINvect(4) - WINvect(3));
  eval([Xname, '=XxXxXx;']);
  eval([Yname, '=YyYyYy;']);
  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);
      for k = perr
	set(sh(k),'string','?');
      end
    end
    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);
      for k = perr
	set(sh(k),'string','?');
      end
    end
    fprintf('\a')
    errordlg(errmsg,'PPLANE error','on');
    return;
  end

  Xname = deblank(sm(1,:));
  Yname = deblank(sm(2,:));

  % If PPLANE Display exists, update it.  If it does not build it.
  
  ppdisp = findobj('name','PPLANE Display');
  if ~isempty(ppdisp)
    figure(ppdisp);
    dh = get(ppdisp,'user');
    ppdispa = dh(18);
    axes(ppdispa);
    cla
    notice = dh(2:3);
    set(notice(1),'string','Computing the field elements.');
    set(notice(2),'string','');
    dh(19:22) = WINvect;
    dh(23:28) = settings;
    set(ppdisp,'user',dh);
    dm = get(dh(1),'user');
    FFFCCNNN = deblank(dm(3,:));
    dm = str2mat(Xname,Yname,FFFCCNNN);
    set(dh(1),'user',dm);
    menu = dh(11:17);
    xmstr = [Xname,' vs. t'];
    ymstr = [Yname,' vs. t'];
    set(menu(3),'label',xmstr);
    set(menu(6),'label',ymstr);
  else
    ppdisp = figure('name','PPLANE Display',...
	'numb','off','visible','off');
    figure(ppdisp);
    Arrflag = settings(1);
    
    % Set up the bulletin window.
    
    nframe = uicontrol('style','frame','units','normal',...
	'pos',[.1 0 .78 .09],'visible','off');
    
    notice(1)=uicontrol('style','text','pos',[.12 .045 .74 .04],...
	'units','normal','horiz','left',...
	'string','Computing the field elements.','visible','off');
    
    notice(2)=uicontrol('style','text','pos',[.12 .005 .74 .04],...
	'units','normal','horiz','left','string','','visible','off');
    
    drawnow
    
    % Set up the buttons and the menu.
    
    global PPRFLAG
    
    RFLAG = 1;
    if size(PPRFLAG) == [1,3]
      RFLAG = 1 - strcmp(PPRFLAG,'off');
    end
    
    
    dbutt(1) = uicontrol('style','push','units','normal','pos',...
	[0.895 0.86 0.09 0.06],'string','Replot',...
	'call','pplane(''refresh'')','visible','off','value',RFLAG);
    
    dbutt(2) = uicontrol('style','push','units','normal','pos',...
	[0.895 0.17 0.09 0.06],'string','Quit',...
	'call','pplane(''quit'')','visible','off');
    
    dbutt(3) = uicontrol('style','push','units','normal','pos',...
	[0.895 0.515 0.09 0.06],'string','Print',...
	'call','pplane(''print'')','visible','off');	
    
    v = version;
    if  (findstr('tudent',v))
      v = 4.2;
    else
      v = str2num(v(1:3));
    end
    dbutt(4) = v;
    
    dbutt(5:7) = [0 0 0];	
    
    menu(1) = uimenu('label','PPLANE Options','visible','off');
    
    menukey = uimenu(menu(1),'label','Keyboard input.','call',...
	'pplane(''kbd'')');
    
    mesev   = uimenu(menu(1),'label','Plot several solutions.',...
	'call','pplane(''several'')');
    
    meplot  = uimenu('label','   Graph','visible','off');
    
    
    menu(3) = uimenu(meplot,'label',[Xname,' vs. t'],...
	'call','pplane(''plotxy'',1)');
    
    menu(6) = uimenu(meplot,'label',[Yname,' vs. t'],...
	'call','pplane(''plotxy'',2)');
    
    meplot3 = uimenu(meplot,'label','Both',...
	'call','pplane(''plotxy'',3)');
    
    meplot4 = uimenu(meplot,'label','3 D',...
	'call','pplane(''plotxy'',4)');
    
    meplot5 = uimenu(meplot,'label','Composite',...
	'call','pplane(''plotxy'',5)');
    
    menu(2) = uimenu(menu(1),'label','Zoom in.',...
	'call','pplane(''zoomin'')','separator','on');
    
    zbmenu = uimenu(menu(1),'label','Zoom back.','call',...
	'pplane(''zoomback'')','enable','off','tag','zbmenu');
    
    
    meeqpt  = uimenu(menu(1),'label','Find an equilibrium point.',...
	'call','pplane(''eqpt'')','separator','on');
    
    menu(4) = uimenu(menu(1),...
	'label','List computed equilibrium points.',...
	'call','pplane(''eqptlist'')');
    
    mestunst= uimenu(menu(1),'label',...
	'Plot stable and unstable orbits.',...
	'call','pplane(''stunst'')');
    
    medall = uimenu(menu(1),'label','Erase all solutions.',...
	'call','pplane(''dall'')','separator','on');
    
    medel = uimenu(menu(1),'label','Delete a graphics object.',...
	'call','pplane(''delete'')','visible','on');
    
    menutext = uimenu(menu(1),...
	'label','Enter text on the Display Window.',...
	'call','pplane(''text'')');
    
    meset = uimenu(menu(1),'label','Settings.',...
	'call','pplane(''settings'')','separator','on');
    
    menu(5) = uimenu(menu(1),...
	'label','Make the Display Window inactive.',...
	'call','pplane(''hotcold'')','separator','on');
    
    if RFLAG
      restring = 'Replot solutions later.';
      buttlist = dbutt([2,3]);
    else
      restring = 'Replot solutions immediately.';
      buttlist = dbutt([1,2,3]);
    end
    
    
    if (v > 4.0)
      menu(7)   = uimenu(menu(1),'label',restring,...
	  'call','pplane(''replot'')','separator','on');
    else
      menu(7) = 0;
    end
    
    if exist('RICE') == 2
      ricemenu = uimenu (menu(1),'label','Set default printer.',...
	  'call','pplane (''riceprint'')','separator','on');
    end
    
    ppdispa = axes('position',[.1 .17 .78 .75],...
	'next','add','box','on',...
	'xgrid','on','ygrid','on','drawmode','fast','visible','off');
    
    set(gcf,'WindowButtonDownFcn','pplane(''down'')');
    
    set(ppdisp,'visible','on');
    set([nframe,notice,buttlist,ppdispa,menu(1),meplot],'visible','on');
    
    % Save the numbers:
    
    % dh(1)			nframe
    % dh(2:3)		notice
    % dh(4:10)		dbutt  (dbutt(5:7) = [0 0 0])
    % dh(11:17)		menu   (menu(7) = 0 if ver = 4.0)
    % dh(18)		ddispa
    % dh(19:22)		WINvect
    % dh(23:28)		settings
    
    dh = [nframe,notice,dbutt,menu,ppdispa,WINvect,settings];
    set(ppdisp,'user',dh);
    
    dm = str2mat(Xname,Yname,'');
    set(dh(1),'user',dm);
    
    printstr = 'print';
    set(dh(11),'user',printstr);
    
  end
  
  % Initialize the window matrix.
  
  set(dh(12),'user',[]);	
  set(findobj('tag','zbmenu'),'enable','off');
  
  
  
  
  % 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
	
	% If an old function m-file exists delete it, and then build a new one.
	
	if (exist(FFFCCNNN)==2) delete([FFFCCNNN,'.m']);end
	tee = clock;
	tee = ceil(tee(6)*100);
	FFFCCNNN=['pptp',num2str(tee)];
	
	fcnstr = ['function YyYypr = ',FFFCCNNN,'(YyYy)\n\n'];
	varstr = [Xname,' = YyYy(1,:);', Yname,' = YyYy(2,:);\n\n'];
	lenstr = ['l = length(YyYy(1,:));\n'];
	derstr1 = ['Xxpr = ', xderivstr,';\n'];
	derstr2 = ['if (length(Xxpr) < l) Xxpr = Xxpr*ones(1,l);end\n'];
	derstr3 = ['Yypr = ', yderivstr,';\n'];
	derstr4 = ['if (length(Yypr) < l) Yypr = Yypr*ones(1,l);end\n'];
	derstr5 = 'YyYypr = [Xxpr;Yypr];\n';
	
	ppf = fopen([FFFCCNNN,'.m'],'w');
	fprintf(ppf,fcnstr);
	fprintf(ppf,varstr);
	fprintf(ppf,lenstr);
	fprintf(ppf,derstr1);
	fprintf(ppf,derstr2);
	fprintf(ppf,derstr3);
	fprintf(ppf,derstr4);
	fprintf(ppf,derstr5);
	fclose(ppf);
	
% Initialize important information as user data.

	dm = str2mat(dm(1:2,:),FFFCCNNN);
	set(dh(1),'user',dm);	% The strings.

	set(dh(2),'user',[]);	% Handles to solution curves 
	                        % and equilibrium points.
	set(dh(3),'user',[]);	% List of equilibrium points.
	set(dh(4),'user',[]);	% Data for solution curves which are 
				% still plotted with no erase mode, 
				% and must be replotted at some point.
	set(dh(5),'user',[]);	% List of equilibrium points which are 
				% still plotted with no erase mode, 
				% and must be replotted at some point.
	set(dh(6),'user',[]);	% Handles for the direction and vector fields.
	set(ppdispa,'user',[]);	% Handles for the ends of orbits plotted with 
				% no erasemode.  These must be deleted at 
				% some point.
	pplane('dirfield');		
	

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 PPLANE Display and get the user data.

	figs = get(0,'children');
	ppdisp = 0;
	for (figno = 1:length(figs))
		if (strcmp(get(figs(figno),'name'),'PPLANE Display'))
			ppdisp = figs(figno);
		end
	end

	dh = get(ppdisp,'user');
	dm = get(dh(1),'user');
	Xname = deblank(dm(1,:));
	Yname = deblank(dm(2,:));
	FFFCCNNN = deblank(dm(3,:));
	
	notice = dh(2:3);
	ppdispa = dh(18);
	WINvect = dh(19:22);
	settings = dh(23:28);

	
	set(notice(1),'string','Computing the field elements.');
	set(notice(2),'string','');
	
% Augment the window matrix

        wmat = get(dh(12),'user');
	wrows = size(wmat,1);
	wflag = 0;
	for k = 1:wrows
	  if wmat(k,:)==WINvect
	    wflag = 1;
	  end
	end
	if wflag == 0
	  wmat = [wmat;WINvect];
	  set(dh(12),'user',wmat);
	end
	if wrows 
	  hhh = findobj('tag','zbmenu');
	  set(hhh,'enable','on');
	end


	
	Xmin = WINvect(1);
	Xmax = WINvect(2);
	Ymin = WINvect(3);
	Ymax = WINvect(4);
	
	N = settings(2);
	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 = feval(FFFCCNNN,[Xx';Yy']);
	Vv = Ww(1,:) + Ww(2,:)*sqrt(-1);
	Vv = Vv.';

	set(notice(1),'string','Computing the graphic elements.');
	set(notice(2),'string','');

	Arrflag = settings(1);
	

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

	arrh = get(dh(6),'user');	% Delete the old field data.
	if (arrh ~= [])
		delete(arrh);
	end

	% We plot both the line field and the vector field.  Then we
	% control which is seen by manipulating the visibility.

	if (Arrflag ==1)
		arrh1 = plot(real(aa1),imag(aa1),'-c','visible','on');
		arrh2 = plot(real(aa2),imag(aa2),'-c','visible','off');
	elseif (Arrflag == 2)
		arrh1 = plot(real(aa1),imag(aa1),'-c','visible','off');
		arrh2 = plot(real(aa2),imag(aa2),'-c','visible','on');
	else
		arrh1 = plot(real(aa1),imag(aa1),'-c','visible','off');
		arrh2 = plot(real(aa2),imag(aa2),'-c','visible','off');
	end
	set(dh(6),'user',[arrh1(:);arrh2(:)]);	% Save the handles for later use.


	% title(Labstr);
	xlabel(Xname);ylabel(Yname);title(' ')
	set(notice(1),'string','Ready.');
	set(notice(2),'string','');
	axis([Dxint,Dyint]);
	pplane('refresh');
	
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);

	dh = 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],...
		'resize','off','name','PPLANE Rice printer selection','numb','off',...
		'vis','off');
	figure(dfrice)
	set(dfrice,'user',dh(11));
	
	
	riceframe = uicontrol('style','frame','pos',frpos,'vis','off');
	printstr = get(dh(11),'user');
	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 = [
	        'dh = get(gcf,''user'');',...
		'me = get(gcf,''currentobject'');',...
		'kk = get(me,''max'');',...
		'set(get(me,''user''),''value'',0),',...
		'set(me,''value'',kk),',...
		'printstr = get(me,''tag'');',...
		'set(dh,''user'',printstr);'];
		
	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(get(gcf,'child'),'units','normal');
	set(gcf,'resize','on');	
	global RICEPRINTPOS
	
	if all(size(RICEPRINTPOS) == [1,4])
	  set(gcf,'pos',RICEPRINTPOS);
	end
	set(gcf,'vis','on');
	set(get(gcf,'child'),'vis','on');





elseif strcmp(action,'replot')

% This is the callback for the Replot button.  This is only visible if
% the version of MATLAB is 4.1 or later.  In this case "enable" is %
% defined and used.  It is only here that RFLAG = get(dbutt(1),'user')
% is changed, so RFLAG will always be equal to 1 if the the version is
% 4.0.

	dh = get(gcf,'user');
	menu = dh(11:17);
	dbutt = dh(4:10);
	RFLAG = strcmp(get(menu(7),'label'), 'Replot solutions immediately.');
	set(dbutt(1),'value',RFLAG);
	if RFLAG
		set(menu(7),'label','Replot solutions later.');
		set(dbutt(1),'visible','off');
	else
		set(menu(7),'label','Replot solutions immediately.');
		set(dbutt(1),'visible','on');
	end
	pplane('refresh')


elseif strcmp(action,'hotcold')

% 'hotcold' is the callback for the menu selection that makes the
% Display Window active or inactive.

	ppdisp = gcf;
	dh = get(ppdisp,'user');
	notice = dh(2:3);
	mehc = dh(15);
	if (findstr(get(mehc,'label'),'inactive'))
		set(ppdisp,'WindowButtonDownFcn',' ');
		set(mehc,'label','Make the Display Window active.');
		set(notice(1),'string','The Display Window is not active.');
		set(notice(2),'string',' ');
	else
		set(ppdisp,'WindowButtonDownFcn','pplane(''down'')');
		set(mehc,'label','Make the Display Window inactive.');
		set(notice(1),'string','Ready.');
		set(notice(2),'string',' ');
	end
	pplane('refresh');
	
elseif strcmp(action,'down')

% 'down' is the Window Button Down call.  It starts the computation of
% solutions from a click of the mouse.

	dh = get(gcf,'user');
	notice = dh(2:3);
	dbutt = dh(4:10);
	RFLAG = get(dbutt(1),'value');
	initpt = get(gca,'currentpoint');
	initpt = initpt(1,[1,2]);
	pplane('solution',initpt);
	if RFLAG
		pplane('refresh');
	end
	set(notice(1),'string','Ready.')
	set(notice(2),'string',' ');
	
	
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.

	dh = get(gcf,'user');
	notice = dh(2:3);
	dbutt = dh(4:10);
	RFLAG = get(dbutt(1),'value');
	set(notice(1),'string','Pick initial points with the mouse.')
	set(notice(2),'string','Enter "Return" when finished.');
	[X,Y]=ginput;
	NN = length(X);
	set(notice(1),'string','Working.')
	set(notice(2),'string',' ');
	for k = 1:NN
		initpt = [X(k),Y(k)];
		pplane('solution',initpt);
	end
	if RFLAG
		pplane('refresh');
	end
	set(notice(1),'string','Ready.')
	set(notice(2),'string',' ');
	
elseif strcmp(action,'solution')

% 'solution' affects the computation and (erasemode == none) plotting of
% solutions.  It also stores the data as appropriate.

	dh = get(gcf,'user');
	dm = get(dh(1),'user');
	notice = dh(2:3);
	RFLAG = get(dh(4),'value');
	if (RFLAG == 0)
		if (strcmp(get(dh(4),'enable'),'off'))
			set(dh(4),'enable','on');
		end
	end
	initpt = input1;
	FFFCCNNN = deblank(dm(3,:));;
	ppdispa = dh(18);
	settings = dh(23:28);
	
	ptstr = [' (',num2str(initpt(1)), ', ', num2str(initpt(2)), ').'];

	set(notice(1),'string',...
		['Working on the forward trajectory from',ptstr])
					
	Dxint = get(ppdispa,'xlim');
	Dyint = get(ppdispa,'ylim');
	htemp = get(ppdispa,'user');
	
	window = [Dxint, Dyint];
	magn = settings(3);
	steps = settings(4);
	itlim = settings(5);
	tol = settings(6);
	
	[tp,xp] = ppsolve(FFFCCNNN,initpt,window,magn,1.0,steps,tol,itlim); 
	hnew = get(ppdispa,'children');
	hnew = hnew(1);
	htemp = [htemp;hnew];	% Save the handles of the orbit ends.
	
	set(notice(1),'string',...
		['Working on the backward trajectory from',ptstr])
	
	[tm,xm] = ppsolve(FFFCCNNN,initpt,window,magn,-1.0,steps,tol,itlim);
	hnew = get(ppdispa,'children');
	hnew = hnew(1);
	htemp = [htemp;hnew];	% Save the handles of the orbit ends.
	
	% Store the trajectory.
	
	xm = flipud(xm);
	x=[xm;xp];
	tm = flipud(tm);
	t=[tm;tp];
	x=[x,t];

	newtraj = get(dh(4),'user');
	[mm,nn] = size(newtraj);
	if (nn == 0)
		newtraj = x;
	else
		% In this case we have to make sure the 
		% trajectories fit in the same matrix.
		l = size(x,1);
		if (l < mm)
			x=[x;ones(mm-l,1)*x(l,:)];
		elseif (mm < l)
			newtraj = [newtraj;ones(l-mm,1)*newtraj(mm,:)];
		end
		newtraj=[newtraj,x];
	end
	set(dh(4),'user',newtraj);
	set(notice(2),'string','Ready.');
	set(ppdispa,'user',htemp);	% Save the handles of the orbit ends.

elseif strcmp(action,'kcompute')

  % 'kcompute' is the call back for the Compute 
  % button on the PPLANE Keyboard figure.

  errmsg1 = 'At least one of the initial values is not a number.';
  errmsg2 = 'You must start PPLANE from the beginning.';
  kh = get(gcf,'user');
  initpt = [str2num(get(kh(2),'string')),str2num(get(kh(4),'string'))];
  if (length(initpt) ~= 2)
    fprintf('\a')
    errordlg(errmsg1,'DFIELD error','on');
    return
  end
  ppdisp = findobj('name','PPLANE Display');
  if isempty(ppdisp)
    fprintf('\a')
    errordlg(errmsg2,'PPLANE error','on');
    return
  else
    figure(ppdisp);
    pplane('solution',initpt);
    dh = get(ppdisp,'user');
    notice = dh(2:3);
    dbutt = dh(4:10);
    RFLAG = get(dbutt(1),'value');
    if RFLAG
      pplane('refresh');
    end
    set(notice(1),'string','Ready.')
    set(notice(2),'string',' ');
  end

		
elseif strcmp(action,'kbd')

  % 'kbd' is the callback for the Keyboard Input menu selection.  It 
  % sets up the PPLANE Keyboard figure which allows accurate input of 
  % initial values using the keyboard.

  dh = get(gcf,'user');
  dm = get(dh(1),'user');
  Xname = deblank(dm(1,:));
  Yname = deblank(dm(2,:));
  
  txtl = 200;txth = 22;
  
  ppkbd = figure('pos',[30 200 (txtl+102) (9*txth)],...
      'resize','off','name','PPLANE Keyboard input','numb','off');
  figure(ppkbd)
  
  KBD(1) = uicontrol('style','text','pos',[1 7*txth txtl txth],...
      'horiz','right','string',['The initial value of ',Xname,' = ']);
  
  KBD(2) = uicontrol('style','edit','pos',[txtl+1 7*txth 100 txth],...
      'string','','call','');
  
  KBD(3) = uicontrol('style','text','pos',[1 5*txth txtl txth],...
      'horiz','right',...
      'string',['The initial value of ',Yname,' = ']);
  
  KBD(4) = uicontrol('style','edit','pos',[txtl+1 5*txth 100 txth],...
      'string','');
  
  KBD(5) = uicontrol('style','push',...
      'pos',[(txtl+102)/2-35 3*txth 70 txth],...
      'string','Compute','call','pplane(''kcompute'')');
  
  KBD(6) = uicontrol('style','push',...
      'pos',[(txtl+102)/2-35,txth,70,txth],...
      'string','Close','call','close');
  
  
  
  set(ppkbd,'user',KBD);
  
  set(ppkbd,'resize','on');
  set(get(ppkbd,'child'),'units','normal');
  
  global PPKBDPOS
  
  if all(size(PPKBDPOS) == [1,4])
    set(gcf,'pos',PPKBDPOS);
  end
  set(gcf,'vis','on');
  set(get(gcf,'child'),'vis','on');
  
	
	
elseif strcmp(action,'eqpt')

% Find and classify equilibrium points.

	dh = get(gcf,'user');
	dm = get(dh(1),'user');
	FFFCCNNN = deblank(dm(3,:));
	notice = dh(2:3);
	dbutt = dh(4:10);
	menu = dh(11:17);
	RFLAG = get(dbutt(1),'value');
	Dx = get(gca,'xlim');
	Dy = get(gca,'ylim');
	DY = [Dx(2)-Dx(1);Dy(2)-Dy(1)];

	epsilon = 1e-4;
	
	set(notice(1),'string','Choose an approximation with the mouse.');
	set(notice(2),'string','');
	z0 = ginput(1);
	
	z = pplane('newton',z0,FFFCCNNN);
	
	flag = z(:,4);
	A = z(:,2:3);
	z = z(:,1);
	
	if (~flag | norm((z-z0')./DY) > 1/5)
		set(notice(1),'string',['There is not an equilibrium point near (',...
			num2str(z0(1)),', ', num2str(z0(2)), ').']);
		set(notice(2),'string','Ready.');
	else
		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

		infostr = str2mat(string1,string2);
		[V,D] = eig(A);			
 		EqPtList = get(dh(3),'user');
		nEqPtList = get(dh(5),'user');

		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];
			set(dh(3),'user',EqPtList);

			if RFLAG
				newh = plot(z(1),z(2),'o','Erasemode','normal');drawnow
				trajh = get(dh(2),'user');
				trajh = [trajh;newh(:)];
				set(dh(2),'user',trajh);
			else
				nEqPtList = [nEqPtList;z',type];
				newh = plot(z(1),z(2),'o','Erasemode','none');drawnow
				htemp = get(gca,'user');
				htemp = [htemp;newh];
				set(gca,'user',htemp);
				set(dh(5),'user',nEqPtList);
				if (strcmp(get(dbutt(1),'enable'),'off'))
					set(dbutt(1),'enable','on');
				end
			end
		end
		pplane('eqptinfo',[V,D,A],infostr);
		set(notice(1),'string','Ready');
		set(notice(2),'string','');
	end

elseif strcmp(action,'eqptinfo')

	dh = get(gcf,'user');
	notice = dh(2:3);
	data = input1;
	V = data(:,1:2);
	D = data(:,3:4);
	A = data(:,5:6);
	infostr = input2; 
	string1 = deblank(infostr(1,:));
	string2 = deblank(infostr(2,:));
	
	ppeqpt = findobj('name','PPLANE Equilibrium point data');
	

	if (ppeqpt)
		eptext =get(ppeqpt,'user');
		set(eptext(1),'string',string1);
		set(eptext(2),'string',string2);
		set(eptext(4),...
			'string',['      ',num2str(A(1,1)),'   ',num2str(A(1,2))]);
		set(eptext(5),...
			'string',['      ',num2str(A(2,1)),'   ',num2str(A(2,2))]);
		set(eptext(7),...
			'string',['     ',num2str(D(1,1)),'    (',num2str(V(1,1)),...
			', ',num2str(V(2,1)),')']);
		set(eptext(8),...
			'string',['     ',num2str(D(2,2)),'    (',num2str(V(1,2)),...
			', ',num2str(V(2,2)),')']);
		figure(ppeqpt)
	else
	
		texth = 20;
		epwind=[0 100 350 12.5*texth];
		epframe=[0 3*texth-3 epwind(3) epwind(4)+3-3*texth];
		ppeqpt = figure('pos',epwind,'resize','off',...
		'name','PPLANE Equilibrium point data','units','pixels',...
		'numb','off');
		figure(ppeqpt)
		epframe = uicontrol('style','frame','pos',epframe);
		goaway = uicontrol('style','push',...
			'pos',[140 texth  70,texth],'string','Go away',...
			'horiz','center','call','close');
			
		eptext(1) = uicontrol('style','text','horiz','left',...
			'pos',[10 11*texth epwind(3)-20 texth],...
			'string',string1);

		eptext(2) = uicontrol('style','text','horiz','left',...
			'pos',[10 10*texth epwind(3)-20 texth],...
			'string',string2);
		
		eptext(3) = uicontrol('style','text','horiz','left',...
			'pos',[10 8.5*texth epwind(3)-20 texth],...
			'string','The Jacobian is:');

		eptext(4) = uicontrol('style','text','horiz','left',...
			'pos',[10 7.5*texth epwind(3)-20 texth],...
			'string',['      ',num2str(A(1,1)),'   ',num2str(A(1,2))]);

		eptext(5) = uicontrol('style','text','horiz','left',...
			'pos',[10 6.5*texth epwind(3)-20 texth],...
			'string',['      ',num2str(A(2,1)),'   ',num2str(A(2,2))]);
			
		eptext(6) = uicontrol('style','text','horiz','left',...
			'pos',[10 5*texth epwind(3)-20 texth],...
			'string','The eigenvalues and eigenvectors are:');

		eptext(7) = uicontrol('style','text','horiz','left',...
			'pos',[10 4*texth epwind(3)-20 texth],...
			'string',['     ',num2str(D(1,1)),'    (',num2str(V(1,1)),...
			', ',num2str(V(2,1)),')']);
		
		eptext(8) = uicontrol('style','text','horiz','left',...
			'pos',[10 3*texth epwind(3)-20 texth],...
			'string',['     ',num2str(D(2,2)),'    (',num2str(V(1,2)),...
			', ',num2str(V(2,2)),')']);
			
		set(ppeqpt,'user',eptext);
		
		set(get(ppeqpt,'child'),'units','normal');
		set(ppeqpt,'resize','on');
	
		global PPEQPTPOS
	
		if all(size(PPEQPTPOS)== [1,4])
		  set(ppeqpt,'pos',PPEQPTPOS);
		end
		set(ppeqpt,'vis','on');
		set(get(ppeqpt,'child'),'vis','on');
		
	end
	
elseif strcmp(action,'stunst')

	dh = get(gcf,'user');
	notice = dh(2:3);
	dbutt = dh(4:10);
	RFLAG = get(dbutt(1),'value');
	Dx = get(gca,'xlim');
	Dy = get(gca,'ylim');
	DY = [Dx(2)-Dx(1);Dy(2)-Dy(1)];
	settings = dh(23:28);
	dm = get(dh(1),'user');
	FFFCCNNN = deblank(dm(3,:));
	EqPtList = get(dh(3),'user');
	nEqPtList = get(dh(5),'user');
	
	Stop = norm(DY)*1e-4;

%  Plot the stable and unstable orbits at a saddle point.

	set(notice(1),'string','Choose a saddle point with the mouse.');
	set(notice(2),'string','');
	z0 = ginput(1);

	z = pplane('newton',z0,FFFCCNNN);
	
	flag = z(:,4);
	A = z(:,2:3);
	z = z(:,1);
	
	if (~flag | norm((z-z0')./DY)> 1/5)
		set(notice(1),'string',['There is not an equilibrium point near (',...
			num2str(z0(1)),', ', num2str(z0(2)), ').']);
	else
		zero = find(abs(z) < 1e-4);
		z(zero) = zeros(size(zero));
		D=det(A);
		if (D>=0)
			set(notice(1),'string',['The equilibrium point at (',...
				num2str(z(1)),', ', num2str(z(2)),...
				') is not a saddle point.']);
		else
			plot(z(1),z(2),'o','Erasemode','none'),drawnow
			htemp = get(gca,'user');
			hnew = get(gca,'children');
			hnew = hnew(1);
			htemp = [htemp;hnew];
			set(gca,'user',htemp);
			if (EqPtList == [])
				EqPtList = [z',1];
				set(dh(3),'user',EqPtList);
				set(dh(5),'user',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];
					nEqPtList = [nEqPtList;z',1];
					set(dh(3),'user',EqPtList);
					set(dh(5),'user',nEqPtList)
				end
			end

			if (RFLAG == 0)
				if strcmp(get(dbutt(1),'enable'),'off')
					set(dbutt(1),'enable','on');
				end
			end

			set(notice(1),'string',['Plotting from the saddle point at (',...
				num2str(z(1)),', ', num2str(z(2)),').']);
			[V,L] = eig(A);
			[L,I] = sort(diag(L));
			
			Veigen = V;
			Ieigen = I;
			Zeigen = z;
			
			window = [Dx, Dy];
			magn = settings(3);
			steps = settings(4);
			itlim = settings(5);
			tol = settings(6);
			offset = norm(DY)/200;	
			lm=abs(L);
			lm=lm/max(lm);
			lm = (max([lm,0.2*ones(size(lm))]'))';
			offs = offset./lm;
			
			
			VeigenV=Veigen(:,Ieigen(1));
			Weigen = Zeigen + offs(1)*VeigenV;
	  	    	set(notice(2),'string','Working on the first stable trajectory.');
		    	[tp,Xeigenp]=ppsolve(FFFCCNNN,Weigen,window,magn,-1,steps,tol,itlim);
			Xeigen1=[Zeigen';Xeigenp];
			x1 = [Xeigen1,[0;tp]];
				
			Weigen = Zeigen - offs(1)*VeigenV;
	  	    	set(notice(2),'string','Working on the second stable trajectory.');
		    	[tp,Xeigenp]=ppsolve(FFFCCNNN,Weigen,window,magn,-1,steps,tol,itlim);
			Xeigen2=[Zeigen';Xeigenp];
			x2 = [Xeigen2,[0;tp]];

			VeigenV=Veigen(:,Ieigen(2));
			Weigen = Zeigen + offs(2)*VeigenV;
	  	    	set(notice(2),'string','Working on the first unstable trajectory.');
		    	[tp,Xeigenp]=ppsolve(FFFCCNNN,Weigen,window,magn,1,steps,tol,itlim);
			Xeigen3=[Zeigen';Xeigenp];
			x3 = [Xeigen3,[0;tp]];
				
			Weigen = Zeigen - offs(2)*VeigenV;
	  	    	set(notice(2),'string','Working on the second unstable trajectory.');
		    	[tp,Xeigenp]=ppsolve(FFFCCNNN,Weigen,window,magn,1,steps,tol,itlim);
			Xeigen4=[Zeigen';Xeigenp];
			x4 = [Xeigen4,[0;tp]];

			l1=size(x1,1);l2=size(x2,1);l3=size(x3,1);l4=size(x4,1);
			mmm = max([l1,l2,l3,l4]);
			x1 = [x1;ones(mmm-l1,1)*x1(l1,:)];
			x2 = [x2;ones(mmm-l2,1)*x2(l2,:)];
			x3 = [x3;ones(mmm-l3,1)*x3(l3,:)];
			x4 = [x4;ones(mmm-l4,1)*x4(l4,:)];
			x  = [x1,x2,x3,x4];

			htemp = get(gca,'user');
			hnew = get(gca,'children');
			hnew = hnew(1:4);
			htemp = [htemp;hnew];
			set(gca,'user',htemp);

			newtraj = get(dh(4),'user');
			[mm,nn] = size(newtraj);
			if (nn == 0)
				newtraj = x;
			else		% In this case we have to make sure the 
					% trajectories fit in the same matrix.
				l = size(x,1);
				if (l < mm)
					x=[x;ones(mm-l,1)*x(l,:)];
				elseif (mm < l)
					newtraj = [newtraj;ones(l-mm,1)*newtraj(mm,:)];
				end
				newtraj=[newtraj,x];
			end
			set(dh(4),'user',newtraj);
			set(ppdispa,'user',htemp);	% Save the handles of the orbit ends.
			if RFLAG
				pplane('refresh');
			end
			set(notice(1),'string','Ready.');
			set(notice(2),'string',' ')
		end
	end

elseif strcmp(action,'zoomin')

% 'zoomin' is the callback for the Zoomin menu item.  It allows the
% user to pick a new display rectangle with two clicks of the mouse.

% This code can not be used until RBBOX works on the Macintosh.
%
%	set(gcf,'WindowButtonDownFcn','pplane(''zoom'')',...
%		'WindowButtonUpFcn','1;','inter','yes');
%	set(gca,'inter','yes');
%	dh = get(gcf,'user');
%	notice = dh(2:3);
%
%	set(notice(1),'string','Pick a new display rectangle by clicking and dragging the mouse.');
%	set(notice(2),'string','');
%
%elseif strcmp(action,'zoom')
%
%	dh = get(gcf,'user');
%	notice = dh(2:3);
%	q1 = get(gcf,'currentpoint');
%	p1 = get(gca,'currentpoint');
%	rbbox([q1 0 0],q1);
%	p2 = get(gca,'currentpoint');
%	a = [p1(1,1:2);p2(1,1:2)];
%	a = [min(a);max(a)];a=a(:);
%	WINvect = a';
%
% This is old code.  To be replaced by that above.
%
	set(gcf,'WindowButtonDownFcn','');
	dh = get(gcf,'user');
	notice = dh(2:3);

	set(notice(1),'string','Pick out two opposite corners with the mouse.');
	set(notice(2),'string','');
	[x1,y1]=ginput(2);
	WINvect = [min(x1), max(x1), min(y1), max(y1)];
%
% This code is good in both cases.
%
	dh(19:22) = WINvect;
	set(gcf,'user',dh);
	figure(gcf);
	axes(gca);
	set(gcf,'WindowButtonDownFcn','pplane(''down'')',...
		'WindowButtonUpFcn','');
	pplane('dirfield');
	ppset = 0;
	figs = get(0,'children');
	for ( figno = 1:length(figs) )
		if(strcmp(get(figs(figno),'name'),'PPLANE Setup'))
			ppset = figs(figno);
		end
	end
	if (ppset)
		sh = get(ppset,'user');
		wind = sh(11:19);
		set(wind(3),'string',num2str(WINvect(1)));
		set(wind(5),'string',num2str(WINvect(2)));
		set(wind(7),'string',num2str(WINvect(3)));
		set(wind(9),'string',num2str(WINvect(4)));
		sh(36:39) = WINvect;
		sh(40:43) = WINvect;
		set(ppset,'user',sh);
	end
	
elseif strcmp(action,'dall')

% 'dall' is the callback for the Erase all solutions option.  

	dh = get(gcf,'user');
	notice = dh(2:3);
	set(notice(1),'string','Erasing all trajectories.');
	set(notice(2),'string','');
	ppdispa = dh(18);
	newtraj = get(dh(4),'user');
	trajh = get(dh(2),'user');
	axes(ppdispa);
	htemp = get(ppdispa,'user');
	if (htemp~=[])
		set(htemp,'erase','normal');
		delete(htemp);			% This deletes the orbit ends.
	end
	set(ppdispa,'user',[]);
	delete(trajh);
	set(dh(2),'user',[]);
	set(dh(3),'user',[]);
	set(dh(4),'user',[]);
	set(dh(5),'user',[]);	
	if (dh(7) > 4.0)
		set(dh(4),'enable','off');
	end
	set(notice(1),'string','Ready.');
	set(notice(2),'string','');
	
elseif strcmp(action,'newton')

% Newton's method to find equilibrium points.

	Iterlimit = 50;
	
	Iter=0;
	
	% The increment for calculating approximate derivatives.
	
	h=.000001;	
	
	FFFCCNNN = input2;
	
	zInit = input1;
	
	zNext=zInit(:);
	functionf=feval(FFFCCNNN,zNext);
	functionf=functionf(:);
	
	% Allow for large/small solutions.
	
	errorlim = norm(functionf,inf)*0.000001;  
	
	while ( (norm(functionf,inf) > errorlim) & (Iter < Iterlimit)	)
		Iter = Iter + 1;
	
	% Now we calculate the jacobian.	
	
		for jjw=1:2			
			sav = zNext(jjw);
			zNext(jjw) = zNext(jjw) + h;
			functionfhj = feval(FFFCCNNN,zNext);	
			functionfhj = functionfhj(:);
			Jacobian(:,jjw) = (functionfhj-functionf)/h;	
			zNext(jjw) = sav;
    	end
    	zNext = zNext - Jacobian\(functionf);
    	functionf = feval(FFFCCNNN,zNext);
		functionf=functionf(:);
	end
	
	if Iter > Iterlimit - 1	    
	    	fLag=[0;0];
	else    		
		fLag=[1;1];	
		for j=1:2			
			sav = zNext(j);
			zNext(j) = zNext(j) + h;
			functionfhj = feval(FFFCCNNN,zNext);	
			functionfhj = functionfhj(:);
			Jacobian(:,j) = (functionfhj-functionf)/h;	
			zNext(j) = sav;
		end
	end
	output = [zNext,Jacobian,fLag];
	
elseif strcmp(action,'paraeval')

%	Replace a parameter with its value in string form.

	para = input1;
	value = input2;
	value = ['(',value,')'];
	str = input3;
	ll = length(str);
	l = length(para);
	lv = length(value);
	if strcmp(para,str)
		str = value;
	elseif (ll >= l+1)
		nstr = [str, 0*para];
		k = 1:ll;
		for j = 1:l
			k = k(find(nstr(k+j-1)==para(j)));
		end
		lk = length(k);
		lopstr = '(+-*/^';
		ropstr = ')+-*/^';
		if (lk>=1)
			s = [];
			pos = 1;
			if ((k(1) == 1) & (find(ropstr == str(l+1))))
				s = value;
				pos = l+1;
			elseif (find(lopstr == str(k(1)-1)))
				if (k(1)+l-1 == ll)
					s =[str(1:(k(1)-1)),value];
					pos = k(1) +l;
				elseif (find(ropstr == str(k(1) + l)))
					s = [str(1:(k(1)-1)),value];
					pos = k(1) +l;
				end
			end
			for j=2:(lk-1)
				if (find(lopstr	== str(k(j)-1)))...
					| (find(ropstr == str(k(j) + l)))
					s = [s,str(pos:(k(j)-1)),value];
					pos = k(j)+l;
				end
			end
			if (lk>1)
				if ( (k(lk) == ll-l+1) & (find(lopstr == str(ll-l))))
					s =[s,str(pos:(k(lk)-1)),value];
					pos = k(lk) + l;
				elseif (find(lopstr	== str(k(lk)-1)))...
						| (find(ropstr == str(k(lk) + l)))
					s =[s,str(pos:(k(lk)-1)),value];
					pos = k(lk) + l;
				end
			end
			str = [s,str(pos:ll)];
		end
	end
	output = str;


elseif strcmp(action,'settings')
	
  % 'settings' is the call back for the Settings menu option.  It sets
  % up the PPLANE Settings window, which allows the user to interactively
  % change several parameters that govern the behaviour of the program.
	
  dh = get(gcf,'user');
  settings = dh(23:28);
  
  txtl = 350;txth = 22;
  figwidth = txtl+104;
  pbwidth = 80;
  framewidth = 3*pbwidth+2*8+2*10;
  frleft = (figwidth - framewidth)/2;
  frbot = 13*txth-5;
  frpos = [frleft frbot framewidth (txth+10)];
  
  ppsettings = figure('pos',[150 50 figwidth (11*txth)],...
      'resize','off','name','PPLANE Settings','numb','off','vis','off');
  figure(ppsettings)
  
  
  setframe = uicontrol('style','frame','pos',frpos);	
  
  settext(1) = uicontrol('style','text',...
      'pos',[3 9*txth txtl txth],'vis','off',...
      'horiz','left','string','Relative size of the calculation window:  ');
		
  kk = 1+2*settings(3);
	
  call2 =[
    'data = get(gcf,''user'');',...
	'nnn = str2num(get(data(7),''string''));',...
	'if isempty(nnn),',...
	'  set(data(7),''string'',''?'');',...
	'  data(3) = NaN;',...
	'else,',...
	'  data(3) = max(0,(nnn-1)/2);',...
	'end,',...
	'set(gcf,''user'',data);'];
  
  setedit(1) = uicontrol('style','edit',...
      'pos',[3+txtl 9*txth 100 txth],...
      'string',num2str(kk),'call',call2,'vis','off');
		
  settext(2) = uicontrol('style','text',...
      'pos',[3 7*txth txtl txth],'vis','off',...
      'horiz','left',...
      'string','Minimum number of plot steps across the window:  ');

  call3 =[
    'data = get(gcf,''user'');',...
	'nnn = str2num(get(data(8),''string''));',...
	'if isempty(nnn),',...
	'  set(data(8),''string'',''?'');',...
	'  data(4) = NaN;',...
	'else,',...
	'  nnn = max(round(nnn),0);',...
	'  set(data(8),''string'',num2str(nnn));',...
	'  data(4) = nnn;',...
	'end,'...
	'set(gcf,''user'',data);'];
		
  setedit(2) = uicontrol('style','edit',...
      'pos',[3+txtl 7*txth 100 txth],...
      'string',num2str(settings(4)),'call',call3,'vis','off');

  settext(3) = uicontrol('style','text',...
      'pos',[3 5*txth txtl txth],'vis','off',...
      'horiz','left','string','Maximum number of iterations per orbit:  ');

  call4 =[
    'data = get(gcf,''user'');',...
	'nnn = str2num(get(data(9),''string''));',...
	'if isempty(nnn),',...
	'  set(data(9),''string'',''?'');',...
	'  data(5) = NaN;',...
	'else,',...
	'  nnn = round(nnn);',...
	'  set(data(9),''string'',num2str(nnn));',...
	'  data(5) = nnn;',...
	'end,'...
	'set(gcf,''user'',data);'];
		
  setedit(3) = uicontrol('style','edit',...
      'pos',[3+txtl 5*txth 100 txth],...
      'string',num2str(settings(5)),...
      'call',call4,'vis','off');

  settext(4) = uicontrol('style','text',...
      'pos',[3 3*txth txtl txth],...
      'horiz','left','string','Relative error tolerance:  ',...
      'vis','off');

  call5 =[
    'data = get(gcf,''user'');',...
	'nnn = str2num(get(data(10),''string''));',...
	'if isempty(nnn),',...
	'  set(data(10),''string'',''?'');',...
	'  data(6) = NaN;',...
	'else,',...
	'  data(6) = nnn;',...
	'  set(data(10),''string'',num2str(nnn));',...
	'end,'...
	'set(gcf,''user'',data);'];
  
  setedit(4) = uicontrol('style','edit',...
      'pos',[3+txtl 3*txth 100 txth],...
      'string',num2str(settings(6)),'call',call5,'vis','off');

  bwidth = 130;
  spacing = (figwidth-2*bwidth-6)/3;
		
  setbutt(1) = uicontrol('style','push',...
      'pos',[3+spacing,txth,bwidth,txth],...
      'string','Cancel','call','close','vis','off');
	
  setbutt(2) = uicontrol('style','push',...
      'pos',[3+2*spacing+bwidth,txth,bwidth,txth],...
      'string','Change settings',...
      'call','pplane(''setchange'')','vis','off');
		
  set(ppsettings,'user',[settings,setedit]);
	
  set(get(ppsettings,'child'),'units','normal');
  set(ppsettings,'vis','on','resize','on');
	
  global PPSETTINGSPOS
	
  if all(size(PPSETTINGSPOS)== [1,4])
    set(ppsettings,'pos',PPSETTINGSPOS);
  end
  set(ppsettings,'vis','on');
  set(get(ppsettings,'child'),'vis','on');


elseif strcmp(action,'setchange')

  % 'setchange' is the callback for the Change button 
  % on the PPLANE Settings window.

  data = get(gcf,'user');
  settings = data(1:6);
  if any(isnan(settings))
    errmsg = 'At least one of the settings is not a number.';
    fprintf('\a')
    errordlg(errmsg,'PPLANE error','on');
    return
  end 
  ppdisp = findobj('name','PPLANE Display');
  if isempty(ppdisp)
    fprintf('\a')
    errordlg('You must start PPLANE from the beginning.',...
	'PPLANE error','on');
  else
    dh = get(ppdisp,'user');
    dh(23:28) = settings;
    set(ppdisp,'user',dh);
  end
  ppset = findobj('name','PPLANE Setup');
  if isempty(ppset)
    fprintf('\a')
    errordlg('You must start PPLANE from the beginning.',...
	'PPLANE error','on');
  else
    sh = get(ppset,'user');
    sh(44:49) = settings;
    set(ppset,'user',sh);
  end	
  close


elseif strcmp(action,'delete')

% 'delete' is the callback for the Delete a graphics object selection
% on the menu.

	oldcall = get(gcf,'WindowButtonDownFcn');
	set(gcf,'WindowButtonDownFcn','');
	dh = get(gcf,'user');
	dm = get(dh(1),'user');
	arrh = get(dh(6),'user');
	v1 = get(arrh(1),'visible');
	v2 = get(arrh(2),'visible');
	set(arrh,'visible','off');
	pplane('refresh')
	trjh = get(dh(2),'user');
	eqpt = get(dh(3),'user');
	notice = dh(2:3);
	set(notice(1),'string','Select a graphics object with the mouse.');
	set(notice(2),'string','');
	ginput(1);
	objh = get(gcf,'currentobject');
	hc = [gcf gca arrh' dh(1:14)] - objh;
	if ( all(hc) ) 
		solh = find(trjh == objh);
		if (solh ~= [])
			trjh = [trjh(1:solh-1);trjh(solh+1:length(trjh))];
			set(dh(2),'user',trjh);
			xd = get(objh,'xdata');
			yd = get(objh,'ydata');
			if length(xd)==1;
				hh = find(eqpt(:,1) == xd);
				if length(hh) > 1
					hh = find(eqpt(hh,2) == yd);
				end
				eqpt = [eqpt(1:hh-1,:);eqpt(hh+1:size(eqpt,1),:)];
				set(dh(3),'user',eqpt);
			end
		end
		delete(objh); 
		set(notice(1),'string','Ready.');
		set(notice(2),'string',' ');
	else
		set(notice(1),'string','The object you selected cannot be deleted.');
		set(notice(2),'string','Ready');
	end
	set(arrh(1),'visible',v1);
	set(arrh(2),'visible',v2);

	set(gcf,'WindowButtonDownFcn',oldcall);

elseif strcmp(action,'refresh')

% 'refresh' is the callback for the Replot the Display Window 
% selection from the menu.  It deletes the ends of orbits which 
% have been plotted using erasemode none, and it replots those 
% solutions so that they become children of the the display 
% window axes.

% 'refresh' is also called at every convenient opportunity in 
% order to keep the Disply Window up to date.

	ppdisp = 0;
	figs = get(0,'children');
	for (figno = 1:length(figs))
		if (strcmp(get(figs(figno),'name'),'PPLANE Display'))
			ppdisp = figs(figno);
		end
	end

	dh = get(ppdisp,'user');
	dm = get(dh(1),'user');
	ppdispa = dh(18);
	newtraj = get(dh(4),'user');
	trajh = get(dh(2),'user');
	axes(ppdispa);
	htemp = get(ppdispa,'user');
	if (htemp~=[])
		set(htemp,'erase','normal');
		delete(htemp);			% This deletes the orbit ends.
	end
	set(ppdispa,'user',[]);
	trjno = size(newtraj,2)/3;
	tees = [1:trjno]*3;
	exes = tees - 2;
	whys = tees - 1;
	newh = line(newtraj(:,exes),newtraj(:,whys),'linestyle','-','color','y');
	for k=1:trjno
		set(newh(k),'Zdata',newtraj(:,3*k));
	end
	EqPtList = get(dh(5),'user');
	if (EqPtList ~= [])
		kk = size(EqPtList,1);
		newhh = zeros(1,kk);
		for k=1:kk
			newhh(k) = plot(EqPtList(k,1),EqPtList(k,2),'o');
		end
	end
	trajh = [trajh;newh(:);newhh(:)];
	drawnow
	set(dh(4),'user',[]);
	set(dh(5),'user',[]);
	set(dh(2),'user',trajh);
	if (dh(7) > 4.0)
		set(dh(4),'enable','off');
	end

elseif strcmp(action,'text')

	dh = get(gcf,'user');
	notice = dh(2:3);
	set(notice(1),'string','Enter the text into the Command Window.');
	set(notice(2),'string',' ');

	oldcall = get(gcf,'WindowButtonDownFcn');
	set(gcf,'WindowButtonDownFcn','');

	txtstr = input('Enter the text here. >  ','s');
	set(notice(1),'string','Choose the location of the text.');
	gtext(txtstr);
	pplane('refresh');
	set(gcf,'WindowButtonDownFcn',oldcall);
	set(notice(1),'string','Ready.');

elseif strcmp(action,'plotxy')

  type = input1;
  ppdisp =gcf;
  ppdispa = gca;
  oldcall = get(ppdisp,'WindowButtonDownFcn');
  set(ppdisp,'WindowButtonDownFcn','');
  xstring = get(get(ppdispa,'XLabel'),'string');
  ystring = get(get(ppdispa,'YLabel'),'string');
  fontn = get(ppdispa,'fontname');
  fonts = get (ppdispa,'fontsize');
  
  
  dh = get(ppdisp,'user');
  dm = get(dh(1),'user');
  notice = dh(2:3);
  handles = get(dh(2),'user');
  sc = get(dh(4),'user');
  printstr = get(dh(11),'user');
  exsol=0;
  if sc~=[]
    exsol=1;
  end
  lll = length(handles);
  
  if (lll >0) & (exsol == 0)
    kk=0;
    while (kk < lll) & (exsol == 0)
      kk = kk + 1;
      xd = get(handles(kk),'xdata');
      if (length(xd) > 7)
	exsol = 1;
      end
    end
  end
  if (exsol == 0)
    set(notice(1),'string','There are no solution curves.');
    set(notice(2),'string','Ready.');
    set(ppdisp,'WindowButtonDownFcn',oldcall);
    return
  end
  
  arrh = get(dh(6),'user');
  v1 = strcmp(get(arrh(1),'visible'),'on');
  v2 = strcmp(get(arrh(2),'visible'),'on');
  set(arrh,'visible','off');
  pplane('refresh')
  set(notice(1),'string','Select a solution curve with the mouse.');
  set(notice(2),'string','');
  ginput(1);
  objh = get(gcf,'currentobject');
  if ( objh == []) 
    solflag = 0;
  elseif ~(strcmp(get(objh,'type'),'line'))
    solflag = 0;
  else
    t=get(objh,'zdata'); 
    if length(t)<10
      solflag = 0;
    else
      solflag = 1;
    end	
  end
  if solflag == 0
    set(notice(1),'string','This object is not a solution curve.');
    set(notice(2),'string','Ready.');
    set(ppdisp,'WindowButtonDownFcn',oldcall);
    if v1 
      set(arrh(1),'visible','on'); 
    elseif v2 
      set(arrh(2),'visible','on'); 
    end
    return
  end
  if v1 
    set(arrh(1),'visible','on'); 
  elseif v2 
    set(arrh(2),'visible','on'); 
  end
  t=get(objh,'zdata');
  x=get(objh,'xdata');
  y=get(objh,'ydata');
  
  texth = 24;
  textw = 100;
  figmargin = 15;
  frmargin = 5;
  sep = 30;
  lwid = textw + 10; 	% The first guess at the width of the legend.
  
  % The absolute position of the t-plot axes
  
  rr = 0.8;
  
  aleft = 73.8; 
  abott = 47.2;
  awidth = 434*rr;
  aht =  342*rr;
  
  frsep = (aht - 10*texth-10)/3;
  
  % The PPLANE t-plot figure --- first guess.
  
  fwid1 = aleft + awidth + sep + figmargin; 
  fwid = fwid1 + lwid;				% This may change.
  fht = abott + aht + 32;	
  
  pos = get(gcf,'position');
  figpos = [pos(1)+30,pos(2)-30,fwid,fht];
  ppxy = figure('name','PPLANE t-plot','position',figpos,'number','off',...
      'visible','off');
  figure(ppxy);
  set(ppxy,'user',[t;x;y]);
  
  
  % The legend.
  
  lleft = aleft + awidth + sep;
  lbot = abott + frsep*3 + texth*8 +10;
  lpos = [lleft/fwid,lbot/fht,lwid/fwid,2*texth/fht];
  
  leg = axes('pos',lpos,'box','on','drawmode','fast','nextplot','add',...
      'xtick',[-1],'ytick',[-1],'xticklabels','','yticklabels','',...
      'xlim',[0,1],'ylim',[0,1],'clipping','on','visible','off');
  
  set(leg,'user',str2mat(xstring,ystring));
  axes(leg);
  xh = text(0,0,xstring,'units','norm','visible','off');
  yh = text(0,0,ystring,'units','norm','visible','off');
  xext = get(xh,'extent');delete(xh);
  yext = get(yh,'extent');delete(yh);
  tw = max(xext(3),yext(3));
  
  % If the legend text is too big, make the legend and the figure larger.
  
  rrrr = 0.6;
  if (tw > rrrr)
    rrr = tw/rrrr;
    lwid = lwid*rrr;
    fwid = fwid1 + lwid;
    figpos(3) = fwid;
    set(ppxy,'pos',figpos);
    lpos(3) = lwid/fwid;
    lpos(1) = lleft/fwid;
    set(leg,'pos',lpos);
    tw = rrrr;
  end
  
  tx = min(1-tw,2/3);
  text('units','norm','pos',[tx,2/3],'string',xstring);
  text('units','norm','pos',[tx,1/3],'string',ystring);
  %	line('xdata',[0,1],'ydata',[0,1],'linestyle','.','visible','off');
  line('xdata',[0.1,tx-0.1],'ydata',[2/3,2/3],'linestyle','-',...
      'color','y');
  line('xdata',[0.1,tx-0.1],'ydata',[1/3,1/3],'linestyle','--',...
      'color','r');	
  
  set(ppxy,'visible','on');
  
  if type ==3 
    set(leg,'visible','on');
    set(get(leg,'children'),'visible','on');
  end
  
  % The t-plot axes.  This changes with fwid.
  
  alft = aleft/fwid; 
  awid = awidth/fwid;
  abot = abott/fht; 
  aht = aht/fht;
  apos = [alft,abot,awid,aht];
  
  axy = axes('pos',apos,'box','on','xgrid','on','ygrid','on','next','add',...
      'drawmode','fast');
  
  % The position of the close button
  
  bbot = abott;
  bwid = lwid;
  bleft = aleft + awidth + sep + (lwid - bwid)/2;
  bpos = [bleft,bbot,bwid,texth];
  
  % The position of the print button.
  
  pbbot = bbot + texth + frsep;
  pbpos = [bleft, pbbot, lwid, texth];
  
  % The positions of the radio frame and its elements.
  
  frbot = pbbot + texth + frsep;
  frleft = lleft;
  frht = 6*texth+10;
  frpos = [frleft,frbot,lwid,frht];
  
  frtitleft = frleft +5;
  frtitwid = lwid -10;
  frtitpos = [frtitleft, frbot+5*texth+5, frtitwid, texth];
  
  radleft = frleft + 5;
  radwid = lwid - 10;
  radbot5 = frbot + 5;
  radbot4 = radbot5 + texth;
  radbot3 = radbot4 + texth;
  radbot2 = radbot3 + texth;
  radbot1 = radbot2 + texth;
  pcall = [
    'ppdisp = findobj(''name'',''PPLANE Display'');',...
	'dh = get(ppdisp,''user'');',...
	'printstr = get(dh(11),''user'');',...
	'eval(printstr)'
    ];
  
  but = uicontrol('style','push','pos',bpos,'string','Go away','call','close');
  
  pbut = uicontrol('style','push','pos',pbpos,'string','Print','call',pcall);
  
  pframe = uicontrol('style','frame','pos',frpos);
  
  pfrtitle = uicontrol('style','text','pos',frtitpos,...
      'string','Graph','horizon','center');
  
  rval1 = ~(type -1);
  rval2 = (~(type -2))*2;
  rval3 = (~(type -3))*3;
  rval4 = (~(type -4))*4;
  rval5 = (~(type -5))*5;
  
  rad(1) = uicontrol('style','radio','pos',[radleft radbot1 radwid texth],...
      'string',xstring,'value',rval1);
  
  rad(2) = uicontrol('style','radio','pos',[radleft radbot2 radwid texth],...
      'string',ystring,'value',rval2,'max',2);
  
  rad(3) = uicontrol('style','radio','pos',[radleft radbot3 radwid texth],...
      'string','Both','value',rval3,'max',3);
  
  rad(4) = uicontrol('style','radio','pos',[radleft radbot4 radwid texth],...
      'string','3 D','value',rval4,'max',4);
  
  rad(5) = uicontrol('style','radio','pos',[radleft radbot5 radwid texth],...
      'string','Composite','value',rval5,'max',5);
  
  for i=1:5
    set(rad(i),'user',[rad(:,[1:(i-1),(i+1):5]),leg,axy]);
  end
  
  callrad = [
    'me = get(gcf,''currentobject'');',...
	'vv = get(me,''value'');'...
	'mm = get(me,''max'');'...
	'if vv ,'...
	'	hand = get(me,''user'');',...
	'	set(hand(1:4),''value'',0);',...
	'	pplane(''plxy'',me);'...
	'end, '...
	'set(me,''value'',mm);'];
  
  
  set(rad,'call',callrad);
  
  set(get(ppxy,'child'),'units','normal');
  set(ppxy,'resize','on');
  
  global PPXYPOS
  
  if all(size(PPXYPOS)== [1,4])
    set(ppxy,'pos',PPXYPOS);
  end
  set(ppxy,'vis','on');
  set(get(ppxy,'child'),'vis','on');
  
  set(notice(1),'string','Ready.');
  set(ppdisp,'WindowButtonDownFcn',oldcall);
  
  pplane('plxy',rad(type))
  
elseif strcmp(action,'plxy')
  
  radbut = input1;
  ud = get(radbut,'user');
  axy = ud(6);
  axis('auto');
  delete(get(axy,'children'));
  leg = ud(5);
  legch = get(leg,'children');
  type = get(radbut,'max');
  if type == 3
    set(leg,'visible','on');
    set(legch,'visible','on');
  else
    set(leg,'visible','off');
    set(legch,'visible','off');
  end
  
  strings = get(leg,'user');
  xstring = deblank(strings(1,:));
  ystring = deblank(strings(2,:));
  data = get(gcf,'user');
  
  t = data(1,:);
  x = data(2,:);
  y = data(3,:);
  
  axes(axy);
  
  if type == 1 
    plot(t,x,'-y');
    view(2);
    xlabel('t')
    ylabel(xstring)
  elseif type == 2
    plot(t,y,'-y');
    view(2);
    xlabel('t')
    ylabel(ystring)
  elseif type == 3
    plot(t,x,'-y',t,y,'--r');
    view(2);
    xlabel('t')
    ylabel([xstring,' and ',ystring])
  elseif type == 4
    plot3(x,y,t);
    view([-30,20]);
    zlabel('t')
    xlabel(xstring);
    ylabel(ystring);
  else
    plot3(x,y,t);
    view([-30,20]);
    zlabel('t')
    xlabel(xstring);
    ylabel(ystring);
    ax = axis;
    plot3(x,y,ax(5)*ones(size(x)),'--');
    plot3(x,ax(4)*ones(size(y)),t,'--');
    plot3(ax(2)*ones(size(x)),y,t,'--');
  end
  
elseif strcmp(action,'print')
  
  dh = get(gcf,'user');
  notice = dh(2:3);
  printstr = get(dh(11),'user');
  set(notice(1),'string','Preparing to print the PPLANE Display Window.');
  set(notice(2),'string','Please be patient.');
  pplane('refresh');
  set(notice(1),'string','Printing the PPLANE Display Window.');
  set(gcf,'pointer','watch');
  eval(printstr);
  set(notice(1),'string','Ready.');
  set(notice(2),'string','');
  set(gcf,'pointer','arrow');
 
elseif strcmp(action, 'eqptlist')

% The labels for equilibrium points.

	disp(' ');
	EqPtType = ['Saddle point.            ';
            	'Nodal sink.              ';
            	'Nodal source.            ';
            	'Spiral sink.             ';
            	'Spiral source.           ';
            	'Spiral equilibrium point.';
            	'Source.                  ';
            	'Sink.                    ';
            	'Unspecified.             ';];
	dh = get(gcf,'user');
	EqPtList = get(dh(3),'user');
	if (EqPtList == [])
		disp('No equilibrium points have been computed.'),disp(' ')
	else
		disp('The following equilibrium points have been calculated:')
		disp(' ')
		L = size(EqPtList,1);
		for k = 1:L
	   		disp([sprintf('(%6.4f, %6.4f)\t',EqPtList(k,1),EqPtList(k,2)),...
			EqPtType(EqPtList(k,3),:)])
		end
		disp(' ')
	end
	
elseif strcmp (action,'zoomback')

  
	dh = get(gcf,'user');
	dm = get(dh(1),'user');
	Xname = deblank(dm(1,:));
	Yname = deblank(dm(2,:));
	wmat = get(dh(12),'user');
	WINvect = dh(19:22);
	
	NN = size(wmat,1);
	txtl = 300;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;
	
	ppzb = figure('pos',[150, 50, figwidth, frbot+frht + 10],...
		'resize','off','name','PPLANE Zoomback','numb','off',...
		'vis','off');
	figure(ppzb)
	set(ppzb,'user',wmat);   	
	
	zbframe = uicontrol('style','frame','pos',frpos,'vis','off');
	wch = 0;j=0;
	while wch == 0
	  j = j+1;
	  if WINvect == wmat(j,:)
	    wch = j;
	  end
	end
	rval = zeros(1,NN);rval(wch) = wch;
	maxval = 1:NN;
	
	winstrs = [];
	for j = 1:NN
	  a = num2str(wmat(j,1));
	  b = num2str(wmat(j,2));
	  c = num2str(wmat(j,3));
	  d = num2str(wmat(j,4));
	  winstrj = [' ',a,' < ',Xname,' < ',b,'   &   ',c,' < ', Yname,' < ',d];
	  if j ==1
	    winstr = winstrj;
	  else
	    winstr = str2mat(winstr,winstrj);
	  end
	end

	for j=1:NN
	  
	  rad(j) = uicontrol('style','radio',...
	      'pos',[(frleft+8), frbot+10+(NN-j)*txth pbwidth+txtl txth],...
	      'string',winstr(j,:),'value',rval(j),'max',j,...
	      'tag',deblank(winstr(j,:)),'hor','left','vis','off');
	end    

			
	for i=1:NN
		set(rad(i),'user',rad(:,[1:(i-1),(i+1):NN]));
	end
	
	set(rad(wch),'tag','me');

	callrad = [
		'me = get(gcf,''currentobject'');',...
		'kk = get(me,''max'');',...
		'set(get(me,''user''),''value'',0),',...
		'set(me,''value'',kk),',...
		'set(get(me,''user''),''tag'','' ''),',...
		'set(me,''tag'',''me''),',...
		];
		
	set(rad,'call',callrad);
	textpos = [frleft+8, frbot+10+(NN+0.5)*txth,pbwidth+txtl,txth];
	textstr = 'Select the rectangle of your choice';
	textw=uicontrol('style','text','pos',textpos,'string',textstr,...
	        'hor','left','vis','off');
	bwidth=80;
	bpos1 = [(figwidth/3-bwidth/2),10,bwidth,txth];
	bpos2 = [(2*figwidth/3-bwidth/2),10,bwidth,txth];
	offbutt = uicontrol('style','push','pos',bpos1,'string','Go away',...
	    'call','delete(gcf)');
	
	zoomcall = [
	    	'wmat = get(gcf,''user'');',...
		'ppset = findobj(''name'',''PPLANE Setup'');',...
	        'ppdisp = findobj(''name'',''PPLANE Display'');',...
		'dh = get(ppdisp,''user'');',...
		'sh = get(ppset,''user'');',...
		'hhh = findobj(gcf,''tag'',''me'');',...
		'kk = get(hhh,''value'');',...
		'delete(gcf),',...
		'WINvect = wmat(kk,:);',...
		'dh(19:22) = WINvect;',...
		'set(ppdisp,''user'',dh);',...
		'wind = sh(11:19);',...
		'set(wind(3),''string'',num2str(WINvect(1)));',...
		'set(wind(5),''string'',num2str(WINvect(2)));',...
		'set(wind(7),''string'',num2str(WINvect(3)));',...
		'set(wind(9),''string'',num2str(WINvect(4)));',...
		'sh(36:39) = WINvect;',...
		'sh(40:43) = WINvect;',...
		'set(ppset,''user'',sh);',...
		'pplane(''dirfield'');'
		];
	    
	procbutt = uicontrol('style','push','pos',bpos2,'string','Zoom',...
	    'call',zoomcall);
	
	
	set(get(gcf,'child'),'units','normal');
	
	global PPZOOMBACKPOS
	
	if all(size(PPZOOMBACKPOS) == [1,4])
	  set(gcf,'pos',PPZOOMBACKPOS);
	end
	set(gcf,'vis','on');
	set(get(gcf,'child'),'vis','on');
	set(gcf,'resize','on');


	
elseif strcmp(action,'quit')

	figs = get(0,'children');
	ppdisp = 0;
	for (figno = 1:length(figs))
		nn = get(figs(figno),'name');
		if (length(nn) > 1)
			if (strcmp(get(figs(figno),'name'),'PPLANE Display'))
				ppdisp = figs(figno);
				dh = get(ppdisp,'user');
				dm = get(dh(1),'user');
				FFFCCNNN = deblank(dm(3,:));
				if (exist(FFFCCNNN)==2) delete([FFFCCNNN,'.m']);end
				close(ppdisp)
			elseif (findstr(get(figs(figno),'name'),'PPLANE')~=[])
				close(figs(figno))		
			end
		end
	end
	if (ppdisp)
	end
end