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

        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(ppset,'user',sethandles);
        hhhh = get(ppset,'child');
        
        % uimenus do not have the 'units' property
        
        kkk = find (hhhh == sysmenu);       
        hhhh = [hhhh(1:kkk-1);hhhh(kkk+1:length(hhhh))];
        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(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.

                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');
        
                        
%       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 (find('*/^' == xderivstr(k+1)))
              xderivstr = [xderivstr(1:k-1), xderivstr(k+1:l)];
            end
          l=l-1;
        end

        l=length(yderivstr);
        for ( k = fliplr(findstr('.',yderivstr)))
            if (find('*/^' == yderivstr(k+1)))
              yderivstr = [yderivstr(1:k-1), yderivstr(k+1:l)];
            end
          l=l-1;
        end
        
%       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','');
        
% 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 = '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);
        
        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'');',...
                '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,'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'');',...
                '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,'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'');',...
                '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,'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'');',...
                '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,'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);
        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');
        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','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;
        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,'--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);
        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
