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: January 31, 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');

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

	pend1 = str2mat(...
		'theta','ang_vel',...
		' ang_vel',' - sin(theta) -D*ang_vel',...
		'D',' ',' ','  ');
	pend2 = str2mat(...
		'0',' ',' ',' ',...
		'2','20','-10','10','-4','4');
	
	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');
	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');
	
	compspec1 = str2mat('x','y',...
		' (1 - x - y)*x',' (A - B*x - C*y)*y',...
		'A','B','C',' ');
	comspec2 = str2mat('4','2','7',' ',...
		'2','20','0','1','0','1');

	vanderpol1 = str2mat('x','y',...
		' M*x - y -x^3',' x',...
		'M',' ',' ',' ');
	vanderpol2 = str2mat('2','0','0','0',...
		'2','20','-5','5','-5','5');

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

	systems1 = str2mat(default1,default2,harmonic1, harmonic2,pend1,pend2,...
		predprey1,predprey2,linear1,linear2); 

	systems2 = str2mat(compspec1,comspec2,vanderpol1,vanderpol2,duffing1,duffing2);

	systems = str2mat(systems1,systems2);

	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:12,:));'...
	'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:12,:));'...
	'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:12,:));'...
	'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:12,:));'...
	'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);'...
	'sh(36)=str2num(get(wind(3),''string''));'...
	'set(gcf,''user'',sh);'];	

	w2 = [
	'sh = get(gcf,''user'');'...
	'wind =sh(11:19);'...
	'sh(37)=str2num(get(wind(5),''string''));'...
	'set(gcf,''user'',sh);'];	

	w3 = [
	'sh = get(gcf,''user'');'...
	'wind =sh(11:19);'...
	'sh(38)=str2num(get(wind(7),''string''));'...
	'set(gcf,''user'',sh);'];	

	w4 = [
	'sh = get(gcf,''user'');'...
	'wind =sh(11:19);'...
	'sh(39)=str2num(get(wind(9),''string''));'...
	'set(gcf,''user'',sh);'];	
	
	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:12,:));'...
	'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:12,:));'...
	'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,:));'...
	'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);'...
	'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);'];
		
	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 = floor(str2num(get(me,''string'')));',...
		'[m,N] = computer;'...
		'if (N <= 8192),',...
		'	N = 32;',...
		'else,',...
		'	N = 50;',...
		'end,'...
		'kk = min([N,max([5,kk])]);'...
		'set(me,''string'',num2str(kk));'...
		'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');

	sysmenu = uimenu('label','Gallery','visible','off');
	
	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(gcf,'user',sethandles);
	set(ppset,'visible','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(frame(1),'user',sm);
	set(frame(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,'system')
	sh = get(gcf,'user');
	set(sh(3),'user',0);
	systems = get(sh(4),'user');
	if strcmp(input1,'default')
		syst = systems(1:18,:);
	elseif strcmp(input1,'harmonic')
		syst = systems(19:36,:);
	elseif strcmp(input1,'pendulum')
		syst = systems(37:54,:);
	elseif strcmp(input1,'prepre')
		syst = systems(55:72,:);
	elseif strcmp(input1,'linear')
		syst = systems(73:90,:);
	elseif strcmp(input1,'compspec')
		syst = systems(91:108,:);
	elseif strcmp(input1,'vanderpol')
		syst = systems(109:126,:);
	elseif strcmp(input1,'duff')
		syst = systems(127:144,:);
	end

	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.

	if (flag == 1)
		WINvect = sh(36:39);
		Arrflag = sh(44);
		NumbFPts = sh(45);
%		set(ppset,'user',sh);	
		ppdisp = 0;
		figs = get(0,'children');
		for (kk = 1:length(figs))
			if (strcmp(get(figs(kk),'name'),'PPLANE Display'))
				ppdisp = figs(kk);
			end
		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;

% Find PPLANE Display if it exists.
	
	figs = get(0,'children');
	ppdisp = 0;
	for (figno = 1:length(figs))
		if (strcmp(get(figs(figno),'name'),'PPLANE Display'))
			ppdisp = figs(figno);
			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','');
		end
	end

% 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);
	
% If PPLANE Display exists, update it.  If it does not build it.
	
	if (ppdisp)
		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.

		dbutt(1) = uicontrol('style','push','units','normal','pos',...
			[0.895 0.86 0.09 0.06],'string','Replot',...
			'call','pplane(''refresh'')','visible','off','value',1);
		
		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');

		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 (v > 4.0)
		 	menu(7)   = uimenu(menu(1),'label','Replot solutions later.',...
			'call','pplane(''replot'')','separator','on');
		else
			menu(7) = 0;
		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,dbutt(2:3),ppdispa,menu(1),meplot],'visible','on');
%		if (dbutt(4) < 4.1)
%			set(menu(7),'visible','off');
%		end

	
	% 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);
		
	end
	
			
% 	Convert the parameters to their current values.
	
%	First remove the blanks.
	xderivstr(find(abs(xderivstr)==32))=[];
	yderivstr(find(abs(yderivstr)==32))=[];
	
%	Replace the parameters with their values.
	if (pname1 ~= [])
		xderivstr = pplane('paraeval',pname1,parav1,xderivstr);
		yderivstr = pplane('paraeval',pname1,parav1,yderivstr);
	end
	if (pname2 ~= [])
		xderivstr = pplane('paraeval',pname2,parav2,xderivstr);
		yderivstr = pplane('paraeval',pname2,parav2,yderivstr);
	end
	if (pname3 ~= [])
		xderivstr = pplane('paraeval',pname3,parav3,xderivstr);
		yderivstr = pplane('paraeval',pname3,parav3,yderivstr);
	end
	if( pname4 ~= [])
		xderivstr = pplane('paraeval',pname4,parav4,xderivstr);
		yderivstr = pplane('paraeval',pname4,parav4,yderivstr);
	end
	
	
% We have to make the derivative strings array smart.
		
	l = length(xderivstr);
	for (k=fliplr(find((xderivstr=='^')|(xderivstr=='*')|(xderivstr=='/'))))
		xderivstr = [xderivstr(1:k-1) '.' xderivstr(k:l)];
		l = l+1;
	end
	l = length(yderivstr);
	for (k=fliplr(find((yderivstr=='^')|(yderivstr=='*')|(yderivstr=='/'))))
		yderivstr = [yderivstr(1:k-1) '.' yderivstr(k:l)];
		l = l+1;
	end
	
% 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','');
	
	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,'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 = 'Enter initial values for both variables.';
	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)
		error(errmsg1);
	end
	figs = get(0,'children');
	ppdisp = 0;
	for (figno = 1:length(figs));
		if (strcmp(get(figs(figno),'name'),'PPLANE Display'))
			ppdisp = figs(figno);
		end
	end
	if (ppdisp == 0)
		error(errmsg2);
	else
		figure(ppdisp);
		pplane('solution',initpt);
	end

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


		
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);
	
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 = 0;
	figs = get(0,'children');
	for (figno = 1:length(figs))
		if (strcmp(get(figs(figno),'name'),'PPLANE Equilibrium point data'))
			ppeqpt = figs(figno);
		end
	end

	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);
	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) == 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');
	figure(ppsettings)
	
	
	setframe = uicontrol('style','frame','pos',frpos);	
	
	settext(1) = uicontrol('style','text','pos',[3 9*txth txtl txth],...
		'horiz','left','string','Relative size of the calculation window:  ');
		
	kk = 1+2*settings(3);
	
	call2 =[
		'data = get(gcf,''user'');',...
		'data(3) = max(0,(str2num(get(data(7),''string''))-1)/2);',...
		'set(gcf,''user'',data);'];
	
		
	setedit(1) = uicontrol('style','edit','pos',[3+txtl 9*txth 100 txth],...
		'string',num2str(kk),'call',call2);
		
	settext(2) = uicontrol('style','text','pos',[3 7*txth txtl txth],...
		'horiz','left','string','Minimum number of plot steps across the window:  ');

	call3 =[
		'data = get(gcf,''user'');',...
		'ss = max(round(str2num(get(data(8),''string''))),0);',...
		'data(4) = ss;',...
		'set(data(8),''string'',num2str(ss));',...
		'set(gcf,''user'',data);'];
		
	setedit(2) = uicontrol('style','edit','pos',[3+txtl 7*txth 100 txth],...
		'string',num2str(settings(4)),'call',call3);

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

	call4 =[
		'data = get(gcf,''user'');',...
		'data(5) = round(str2num(get(data(9),''string'')));',...
		'set(data(9),''string'',num2str(data(5)));',...
		'set(gcf,''user'',data);'];
		
	setedit(3) = uicontrol('style','edit','pos',[3+txtl 5*txth 100 txth],...
		'string',num2str(settings(5)),'call',call4);

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

	call5 =[
		'data = get(gcf,''user'');',...
		'data(6) = str2num(get(data(10),''string''));',...
		'set(data(10),''string'',num2str(data(6)));',...
		'set(gcf,''user'',data);'];
		
	setedit(4) = uicontrol('style','edit','pos',[3+txtl 3*txth 100 txth],...
		'string',num2str(settings(6)),'call',call5);

	bwidth = 130;
	spacing = (figwidth-2*bwidth-6)/3;
		
	setbutt(1) = uicontrol('style','push','pos',[3+spacing,txth,bwidth,txth],...
		'string','Cancel','call','close');
	
	setbutt(2) = uicontrol('style','push','pos',[3+2*spacing+bwidth,txth,bwidth,txth],...
		'string','Change settings','call','pplane(''setchange'')');
		
	set(ppsettings,'user',[settings,setedit]);

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);
	ppdisp = 0;
	figs = get(0,'children');
	for (figno = 1:length(figs));
		if (strcmp(get(figs(figno),'name'),'PPLANE Display')),
			ppdisp = figs(figno);
		end
	end
	if (ppdisp == 0)
		error('You must start PPLANE from the beginning.');
	else
		dh = get(ppdisp,'user');
		dh(23:28) = settings;
		set(ppdisp,'user',dh);
	end
	ppset = 0;
	for (figno = 1:length(figs));
		if (strcmp(get(figs(figno),'name'),'PPLANE Setup'))
			ppset = figs(figno);
		end
	end
	if (ppset == 0)
		error('You must start PPLANE from the beginning.');
	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');
	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','b');	

	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;

	but = uicontrol('style','push','pos',bpos,'string','Go away','call','close');

	pbut = uicontrol('style','push','pos',pbpos,'string','Print','call','print');

	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(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,'--b');
		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
%	grid


elseif strcmp(action,'print')

	dh = get(gcf,'user');
	notice = dh(2:3);
	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');
	print 
	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,'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
