matpplane.m

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

%   matpplane  is an interactive tool for studying planar autonomous
%	systems of differential equations.  When matpplane is executed, a
%	matpplane Setupwindow 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 matpplane
%	Display window is opened, and a field of the type requested is
%	displayed for the system.  When the mouse button is depressed in
%	the matpplane Display window, the solution to the system with that
%	initial condition is calculated and plotted.
%
%	Other options are available in the menus.  These are
%	fairly self-explanatory.
%


%  Copyright (c)  1995, 1997, 1998, 1999, 2000, 2001, 2002, 2003
%                  John C. Polking, Rice University
%                  Last Modified: April 24, 2003

%  Modified by JDM Nov 2016 for Matlab 2016.
%  Modified by Nathanael Kazmierczak Jan 2021 for Matlab 2021.
%  Modified by Roy Goodman Jan 2021 for Matlab 2020b


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



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

if strcmp(action,'initialize')
    
    % First thing that runs when the program starts
    % First we make sure that there is no other copy of matpplane
    % running, since this causes problems.
    pph = findobj('tag','matpplane');
    if ~isempty(pph)
        qstring = {'There are some matpplane figures open although they may be invisible.  ';...
            'What do you want to do?'};
        tstring = 'Only one copy of matpplane can be open at one time.';
        b1str = 'Restart matpplane.';
        b2str = 'Just close those figures.';
        b3str = 'Do nothing.';
        answer = questdlg(qstring,tstring,b1str,b2str,b3str,b1str);
        if strcmp(answer,b1str)
            delete(pph);
            %matpplane;return
            eval(startstr);return
        elseif strcmp(answer,b2str)
            delete(pph);return
        else
            return
        end
    end  % if ~isempty(pph);
    
    % Make sure tempdir is on the MATLABPATH.  We want to be sure that we
    % return the path to its current position when we exit.
    

%    p = path;
%    tmpdir = tempdir;
%    ll = length(tmpdir);
%    tmpdir = tmpdir(1:ll-1);
    ud.remtd = 0;
    % The following is a more robust way to check if the operating system's
    % temporary directory is on the path
    pathCell = regexp(path, pathsep, 'split');
    if ispc  % Windows is not case-sensitive
        onPath = any(strcmpi(tempdir, pathCell));
    else
        onPath = any(strcmp(tempdir, pathCell));
    end
    
    if ~onPath
        ud.remtd = 1;
        addpath(tempdir)
    end
    
    % Next we look for old files created by matpplane.
    % First in the current directory.
    
    oldfiles = dir('pptp*.m');
    ofiles = cell(0,1);
    % kk = zeros(0,1);
    for k = 1:length(oldfiles)
        fn = oldfiles(k).name;
        fid = fopen(fn,'r');
        fgetl(fid);   %read and discard lines
        fgetl(fid);
        ll = fgetl(fid);
        fclose(fid);
        if strcmp(ll,'%% Created by matpplane')
            delete(fn)
        else
            l= length(ofiles);
            ofiles{l+1} = fn;
        end
    end
    
    %Then in the temp directory.
    
    oldfiles = dir([tempdir,'pptp*.m']);
    for k = 1:length(oldfiles)
        fn = [tempdir,oldfiles(k).name];
        fid = fopen(fn,'r');
        fgetl(fid);    %read and discard lines
        fgetl(fid);
        ll = fgetl(fid);
        fclose(fid);
        if strcmp(ll,'%% Created by matpplane')
            delete(fn)
        else
            l= length(ofiles);
            ofiles{l+1} = fn;
        end
    end
    lll = length(ofiles);
    if lll >0
        if lll == 1
            astr = 'The file';
            bstr = 'If so it can be safely deleted.';
        else
            astr = 'The files';
            bstr = 'If so they can be safely deleted.';
        end
        fprintf([astr,'\n']);
        for j = 1:lll
            fn = ofiles{j};
            disp(['     ',fn]);
        end
        fprintf('may have been created by an older version of PPLANE.\n');
        fprintf(bstr,'\n\n');
    end
    
    style = 'white';
    ppdir = pwd;
    ssize = get(0,'defaultaxesfontsize');
    npts = 20;
    solver = 'Dormand Prince';
    tol = 1e-4;
    stepsize = 0.1;
    if exist('ppstart','file')
        H = ppstart;
        if ~isempty(H)
            if isfield(H,'style')
                style = H.style;
            end
            if isfield(H,'size')
                ssize = H.size;
            end
            if isfield(H,'npts')
                npts = H.npts;
            end
            if isfield(H,'solver')
                solver = H.solver;
            end
            if isfield(H,'ppdir')
                ppdir = H.ppdir;
            end
            if isfield(H,'stepsize')
                stepsize = H.stepsize;
            end
            if isfield(H,'tolerance')
                tol = H.tolerance;
            end
        end
    end
    if get(0,'screendepth') == 1
        style = 'bw';
    end
    
    ud.ssize = ssize;
    ud.ppdir = ppdir;
    comp = computer;
    if strcmp(comp,'PCWIN')
        ud.fontsize = 0.8*ud.ssize;
    else
        ud.fontsize = ud.ssize;
    end
    
    % Set up for the menu of systems.
    
    system.name = 'default system';
    system.xvar = 'x';
    system.yvar = 'y';
    system.xder = ' 2*x - y + 3*(x^2-y^2) + 2*x*y';
    system.yder = ' x - 3*y - 3*(x^2-y^2) + 3*x*y';
    system.pname = {};
    system.pval = {};
    system.fieldtype = 'arrows';
    system.npts = npts;
    system.wind = [-2 4 -4 2];
    
    system(2).name = 'linear system';
    system(2).xvar = 'x';
    system(2).yvar = 'y';
    system(2).xder = ' A*x + B*y';
    system(2).yder = ' C*x + D*y';
    system(2).pname = {'A', 'B', 'C', 'D', '', ''};
    system(2).pval = {'2', '2', '-2',  '-3'};
    system(2).fieldtype = 'arrows';
    system(2).npts = npts;
    system(2).wind = [-5 5 -5 5];
    
    system(3).name = 'vibrating spring';
    system(3).xvar = 'x';
    system(3).yvar = 'v';
    system(3).xder = ' v';
    system(3).yder = ' -(k*x + d*v)/m';
    system(3).pname = {'k', 'm', 'd', '', '', ''};
    system(3).pval = {'3', '1', '0'};
    system(3).fieldtype = 'arrows';
    system(3).npts = npts;
    system(3).wind = [-5 5 -5 5];
    
    system(4).name = 'pendulum';
    system(4).xvar = '\theta';
    system(4).yvar = '\omega';
    system(4).xder = ' \omega';
    system(4).yder = ' -sin(\theta) - D*\omega';
    system(4).pname = {'D', '', '', '', '', ''};
    system(4).pval = {'0'};
    system(4).fieldtype = 'arrows';
    system(4).npts = npts;
    system(4).wind = [-10 10 -4 4];
    
    system(5).name = 'predator prey';
    system(5).xvar = 'prey';
    system(5).yvar = 'predator';
    system(5).xder = ' (A - B*predator)*prey';
    system(5).yder = ' (D*prey - C)*predator';
    system(5).pname = {'A', 'B', 'C', 'D', '', ''};
    system(5).pval = {'0.4', '0.01', '0.3', '0.005'};
    system(5).fieldtype = 'arrows';
    system(5).npts = npts;
    system(5).wind = [0 120 0 80];
    
    system(6).name = 'competing species';
    system(6).xvar = 'x';
    system(6).yvar = 'y';
    system(6).xder = ' r*(1 - x - A*y)*x';
    system(6).yder = ' s*(1 - y - B*x)*y';
    system(6).pname = {'r', 's', 'A', 'B', '', ''};
    system(6).pval = {'0.4', '0.6', '5', '4'};
    system(6).fieldtype = 'nullclines';
    system(6).npts = npts;
    system(6).wind = [0 1 0 1];
    
    system(7).name = 'cooperative species';
    system(7).xvar = 'x';
    system(7).yvar = 'y';
    system(7).xder = ' r*(1 - x + A*y)*x';
    system(7).yder = ' s*(1 - y + B*x)*y';
    system(7).pname = {'r', 's', 'A', 'B', '', ''};
    system(7).pval = {'0.4', '0.6', '2', '0.3'};
    system(7).fieldtype = 'nullclines';
    system(7).npts = npts;
    system(7).wind = [0 9 0 5];
    
    system(8).name = 'van der Pol''s equation';
    system(8).xvar = 'x';
    system(8).yvar = 'y';
    system(8).xder = ' M*x - y - x^3';
    system(8).yder = ' x';
    system(8).pname = {'M', '', '', '', '', ''};
    system(8).pval = {'2'};
    system(8).fieldtype = 'arrows';
    system(8).npts = npts;
    system(8).wind = [-5 5 -5 5];
    
    system(9).name = 'Duffing''s equation';
    system(9).xvar = 'x';
    system(9).yvar = 'y';
    system(9).xder = ' y';
    system(9).yder = ' -(k*x + c*y + l*x^3)/m';
    system(9).pname = {'k', 'c', 'm', 'l', '', ''};
    system(9).pval = {'-1', '0', '1', '1'};
    system(9).fieldtype = 'arrows';
    system(9).npts = npts;
    system(9).wind = [-3 3 -3 3];
    
    system(10).name = 'square limit set';
    system(10).xvar = 'x';
    system(10).yvar = 'y';
    system(10).xder = ' (y + x/5)*(1-x^2)';
    system(10).yder = ' -x*(1-y^2)';
    system(10).pname = {'', '', '', '', '', ''};
    system(10).pval = {'', '', '', ''};
    system(10).fieldtype = 'arrows';
    system(10).npts = npts;
    system(10).wind = [-1.5 1.5 -1 1];
    
    
    ud.c = system(1);	% Changed values.
    pname = ud.c.pname;
    for kk = length(pname)+1:6
        pname{kk} = '';
    end
    pval = ud.c.pval;
    for kk = length(pval)+1:6
        pval{kk} = '';
    end
    ud.c.pname = pname;
    ud.c.pval = pval;
    ud.o = ud.c;	% Original values.
    % ud.h = system(1);	% This will be the handles in the
    % setup window.
    
    ud.style = style;
    ud.npts = npts;
    ud.settings.magn = 1.25;
    ud.settings.refine = 8;
    ud.settings.tol = tol;
    ud.settings.solver = solver;
    ud.settings.stepsize = stepsize;
    ud.settings.speed = 100;
    ud.system = system;
    ud.solvers = {'ppdp45';
        'rk4';
        'ode45';
        'ode23';
        'ode113';
        'ode15s';
        'ode23s';
        'ode23t';
        'ode23tb'};
    
    switch style  % Note that an earlier line explicitly sets style = 'white' so none of this gets called
        case 'black'
            color.temp = [1 0 0]; % red for temporary orbits
            color.orb = [1 1 0];  % yellow for orbits
            color.eqpt = [1 0 0];  % red for eq. pts.
            color.arrows = .5*[1 1 1];  % gray for arrows
            color.narrows = .7*[1 1 1];  % gray for nullcline arrows
            color.tx = [1 1 0]; % yellow for xt plots & 3D plots
            color.ty = [1 0 0]; % red for yt plots
            color.ta = [1 0 0]; % red for axis plots
            color.sep = [.2,1,0]; % green for separatrices
            color.xcline = [1 1 1]; % white for xclines
            color.ycline = [1 0 1]; % magenta for yclines
            color.level = [1,.5,.5];
            
        case 'white'
            color.temp = [1 0 0]; % red for temporary orbits
            color.orb = [0 0 1];  % blue for orbits
            color.eqpt = [1 0 0];  % red for eq. pts.
            color.arrows = 0.7*[1 1 1]; % gray for arrows
            color.narrows = .4*[1 1 1];  % gray for nullcline arrows
            color.tx = [0 0 1]; % blue for xt plots & 3D plots
            color.ty = [1 0 0]; % red for yt plots
            color.ta = [1 0 0]; % red for axis plots
            color.sep = [0,1,0];% green for separatrices
            color.xcline = [1 0 .75]; % magenta for xclines
            color.ycline = [1 .5 0]; %  orange for yclines
            color.level = 0.8*[.9,.5,.8];
            
        case 'test'
            color.temp = [1 0 0]; % red for temporary orbits
            color.orb = [0 0 1];  % blue for orbits
            color.arrows = .7*[1 1 1]; % gray for arrows
            color.eqpt = [1 0 0];  % red for eq. pts.
            color.narrows = .4*[1 1 1];  % gray for nullcline arrows
            color.tx = [0 0 1]; % blue for xt plots & 3D plots
            color.ty = [1 0 0]; % red for yt plots
            color.ta = [1 0 0]; % red for axis plots
            color.sep = [0,1,0];% green for separatrices
            color.xcline = [1 0 .75]; % magenta for xclines
            color.ycline = [1 .5 0]; %  orange for yclines
            color.level = [1,.5,.5];
            
        case 'display'
            color.temp = [1 0 0]; % red for temporary orbits
            color.orb = [0 0 1];  % blue for orbits
            color.arrows = .4*[1 1 1]; % gray for arrows
            color.eqpt = [1 0 0];  % red for eq. pts.
            color.narrows = .4*[1 1 1];  % gray for nullcline arrows
            color.tx = [0 0 1]; % blue for xt plots & 3D plots
            color.ty = [1 0 0]; % red for yt plots
            color.ta = [1 0 0]; % red for axis plots
            color.sep = 2*[.5 0 .5];% [0,1,0];% green for separatrices
            color.xcline = [1 0 .75]; % magenta for xclines
            color.ycline = [1 .5 0]; %  orange for yclines
            color.level = 0.8*[.9,.5,.8];
            
        case 'bw'
            color.temp = [1 1 1]; % white for everything
            color.orb = [1 1 1];
            color.eqpt = [1 1 1];
            color.arrows = [1 1 1];
            color.narrows = [1 1 1];
            color.tx = [1 1 1];
            color.ty = [1 1 1];
            color.ta = [1 1 1];
            color.sep = [1 1 1];
            color.xcline = [1 1 1];
            color.ycline = [1 1 1];
            color.level = [1,1,1];
    end
    ud.color = color;
    ud.level = ' ';
    ppset = figure('name','matpplane Setup',...
        'tag','matpplane','visible','off',...
        'UserData',ud);
    %  ppset = figure('name','matpplane Setup','numb','off',...
    %		 'tag','matpplane','visible','off',...
    %		 'UserData',ud);
    
    matpplane('figdefault',ppset);
    frame(1) = uicontrol('style','frame','visible','off');
    eq(1)=uicontrol('style','text',...
        'horizon','center',...
        'string','The differential equations.','visible','off');
    ext = get(eq(1),'extent');
    rr=ext(4)/10;
    
    texth =ext(4)+4;    % 19;   		% Height of text boxes.
    varw = 40*rr;		% Length of variable boxes.
    equalw =13*rr;		% Length of equals.(30)
    eqlength = 230*rr;		% Length of right hand sides of equations.
    winstrlen = 120*rr;	% Length of string boxes in display frame.
    left = 0;		% Left margin of the frames.
    frsep = 1;    %3;		% Separation between frames.
    separation = texth;	% Separation between boxes.
    
    dfigwidth =2*left + varw+equalw+eqlength+10;	% Width of the figure.
    dfigurebot = 30;	% Bottom of the figure.
    buttw = dfigwidth/3;
    qwind = [0,frsep,buttw,separation];	% Quit button
    rwind = [buttw,frsep,buttw,separation];	% Revert "
    pwind = [2*buttw,frsep,buttw,separation];	% Proceed "
    
    disfrbot = 2*frsep + separation;	% Display frame.
    disfrw = winstrlen + varw +10;
    disfrht = 5*separation + 10;
    disfrwind = [left, disfrbot, disfrw, disfrht];
    
    pfrbot = disfrbot + disfrht +frsep;	% Parameter frame.
    pfrw = dfigwidth -2*left;
    pfrht = 3*separation + 10;
    pfrwind = [left, pfrbot, pfrw, pfrht];
    
    defrbot = pfrbot + pfrht + frsep;	% Equation frame.
    defrw = pfrw;
    defrht = 3*separation + 10;
    defrwind = [left, defrbot, defrw, defrht];
    
    ffrbot = disfrbot;			% Field frame.
    ffrleft = left + disfrw + frsep;
    ffrw = dfigwidth -left - ffrleft;
    ffrht = disfrht;
    ffrwind = [ffrleft, ffrbot, ffrw, ffrht];
    
    dfigureheight = defrbot + defrht +frsep;	% Height of the figure.
    
    set(ppset,'pos',[30 dfigurebot dfigwidth dfigureheight]);
    
    set(frame(1),'pos',defrwind);
    
    
    xname=[
        'ud = get(gcf,''UserData'');'...
        'Xname=get(ud.h.xvar,''string'');'...
        'minxstr = [''The minimum value of '',Xname,'' = ''];',...
        'set(ud.h.twind(1),''string'',minxstr);'...
        'maxxstr = [''The maximum value of '',Xname,'' = ''];',...
        'set(ud.h.twind(2),''string'',maxxstr);'...
        'ud.c.xvar = Xname;'...
        'ud.flag = 0;'...
        'ud.c.name = '''';',...
        'set(gcf,''UserData'',ud);'];
    
    yname=[
        'ud = get(gcf,''UserData'');'...
        'Yname=get(ud.h.yvar,''string'');'...
        'minystr = [''The minimum value of '',Yname,'' = ''];',...
        'set(ud.h.twind(3),''string'',minystr);'...
        'maxystr = [''The maximum value of '',Yname,'' = ''];',...
        'set(ud.h.twind(4),''string'',maxystr);'...
        'ud.c.yvar = Yname;'...
        'ud.flag = 0;'...
        'ud.c.name = '''';',...
        'set(gcf,''UserData'',ud);'];
    
    xder =[
        'ud = get(gcf,''UserData'');'...
        'ud.c.xder = get(ud.h.xder,''string'');'...
        'ud.flag = 0;'...
        'ud.c.name = '''';',...
        'set(gcf,''UserData'',ud);'];
    
    yder =[
        'ud = get(gcf,''UserData'');'...
        'ud.c.yder = get(ud.h.yder,''string'');'...
        'ud.flag = 0;'...
        'ud.c.name = '''';',...
        'set(gcf,''UserData'',ud);'];
    
    equationbot = defrbot + 5;
    eqlabelbot = equationbot + 2*separation;
    xbot = equationbot + separation;		% Bottom of x equation.
    ybot = equationbot;				% Bottom of y equation.
    
    lablen =200*rr;
    
    eqlableft = (dfigwidth-lablen)/2;
    eqleft = left + 5;
    fudge = 0.15*separation;
    
    set(eq(1),'pos',[eqlableft eqlabelbot lablen texth]);
    tcolor = get(gcf,'defaultuicontrolbackgroundcolor');
    ecolor = 'w';
    
    ud.h.xvar=uicontrol('pos',[eqleft, xbot, varw, texth],...
        'style','edit',...
        'horizon','right',...
        'string',ud.o.xvar,...
        'CallBack',xname,...
        'backgroundcolor',ecolor,...
        'visible','off');
    
    [~] = uicontrol('style','text',...
        'pos',[eqleft+varw xbot equalw texth],...
        'horizon','center',...
        'string',''' = ',...
        'backgroundcolor',tcolor,...
        'visible','off');
    
    ud.h.xder=uicontrol('pos',[eqleft+varw+equalw xbot eqlength texth],...
        'string',ud.o.xder,...
        'horizon','left','style','edit',...
        'backgroundcolor',ecolor,...
        'CallBack',xder,'visible','off');
    
    ud.h.yvar=uicontrol('pos',[eqleft ybot varw texth],...
        'style','edit',...
        'horizon','right',...
        'string',ud.o.yvar,...
        'backgroundcolor',ecolor,...
        'CallBack',yname,'visible','off');
    
    [~] = uicontrol('style','text',...
        'pos',[eqleft+varw ybot equalw texth],...
        'horizon','center','string',''' = ',...
        'backgroundcolor',tcolor',...
        'visible','off');
    
    ud.h.yder=uicontrol('pos',[eqleft+varw + equalw ybot eqlength texth],...
        'string',ud.o.yder,...
        'horizon','left','style','edit',...
        'backgroundcolor',ecolor,...
        'CallBack',yder,'visible','off');
    
    
    [~] = uicontrol('style','frame','pos',disfrwind,'visible','off');
    
    w1 = [
        'ud = get(gcf,''UserData'');'...
        'nnn = 	str2num(get(ud.h.wind(1),''string''));'...
        'if isempty(nnn),',...
        '  set(ud.h.wind(1),''string'',''?'');',...
        '  nnn = NaN;',...
        'end,',...
        'ud.c.wind(1) = nnn;',...
        'ud.c.name = '''';',...
        'set(gcf,''UserData'',ud);'];
    
    w2 = [
        'ud = get(gcf,''UserData'');'...
        'nnn = 	str2num(get(ud.h.wind(2),''string''));'...
        'if isempty(nnn),',...
        '  set(ud.h.wind(2),''string'',''?'');',...
        '  nnn = NaN;',...
        'end,',...
        'ud.c.wind(2) = nnn;',...
        'ud.c.name = '''';',...
        'set(gcf,''UserData'',ud);'];
    
    w3 = [
        'ud = get(gcf,''UserData'');'...
        'nnn = 	str2num(get(ud.h.wind(3),''string''));'...
        'if isempty(nnn),',...
        '  set(ud.h.wind(3),''string'',''?'');',...
        '  nnn = NaN;',...
        'end,',...
        'ud.c.wind(3) = nnn;',...
        'ud.c.name = '''';',...
        'set(gcf,''UserData'',ud);'];
    
    w4 = [
        'ud = get(gcf,''UserData'');'...
        'nnn = 	str2num(get(ud.h.wind(4),''string''));'...
        'if isempty(nnn),',...
        '  set(ud.h.wind(4),''string'',''?'');',...
        '  nnn = NaN;',...
        'end,',...
        'ud.c.wind(4) = nnn;',...
        'ud.c.name = '''';',...
        'set(gcf,''UserData'',ud);'];
    
    
    winbot1 = disfrbot + disfrht - 5 - separation;
    winbot2 = winbot1 - separation;
    winbot3 = winbot2 - separation;
    winbot4 = winbot3 - separation;
    winbot5 = winbot4 - separation;
    
    [~] = uicontrol('style','text',...
        'pos',[eqleft winbot1 disfrw-10 texth],...
        'horizon','center',...
        'string','The display window.','visible','off');
    
    % ud.h.twind contains the handles to the text windows, and ud.h.wind
    % contains the handles to the edit windows.
    
    twstr1 = ['The minimum value of ',ud.o.xvar,' = '];
    ud.h.twind(1) = uicontrol('style','text',...
        'pos',[eqleft winbot2 winstrlen texth],...
        'horizon','right',...
        'string',twstr1,...
        'backgroundcolor',tcolor',...
        'visible','off');
    
    ud.h.wind(1) = uicontrol('style','edit',...
        'pos',[eqleft+winstrlen winbot2 40*rr texth],...
        'string',num2str(ud.o.wind(1)),...
        'backgroundcolor',ecolor,...
        'CallBack',w1,'visible','off');
    
    twstr2 = ['The maximum value of ',ud.o.xvar,' = '];
    ud.h.twind(2) = uicontrol('style','text',...
        'pos',[eqleft winbot3 winstrlen texth],...
        'horizon','right',...
        'string',twstr2,...
        'backgroundcolor',tcolor',...
        'visible','off');
    
    ud.h.wind(2) = uicontrol('style','edit',...
        'pos',[eqleft+winstrlen winbot3 40*rr texth],...
        'string',num2str(ud.o.wind(2)),...
        'backgroundcolor',ecolor,...
        'CallBack',w2,'visible','off');
    
    twstr3 = ['The minimum value of ',ud.o.yvar,' = '];
    ud.h.twind(3)= uicontrol('style','text',...
        'pos',[eqleft winbot4 winstrlen texth],...
        'horizon','right',...
        'string',twstr3,...
        'backgroundcolor',tcolor',...
        'visible','off');
    
    ud.h.wind(3) = uicontrol('style','edit',...
        'pos',[eqleft+winstrlen winbot4 40*rr texth],...
        'string',num2str(ud.o.wind(3)),...
        'backgroundcolor',ecolor,...
        'CallBack',w3,'visible','off');
    
    twstr4 = ['The maximum value of ',ud.o.yvar,' = '];
    ud.h.twind(4) = uicontrol('style','text',...
        'pos',[eqleft winbot5 winstrlen texth],...
        'horizon','right',...
        'string',twstr4,...
        'backgroundcolor',tcolor',...
        'visible','off');
    
    ud.h.wind(4) = uicontrol('style','edit',...
        'pos',[eqleft+winstrlen winbot5 40*rr texth],...
        'string',num2str(ud.o.wind(4)),...
        'backgroundcolor',ecolor,...
        'CallBack',w4,'visible','off');
    
    
    [~]=uicontrol('style','frame','pos',pfrwind,'visible','off');
    
    pnCallBack = [
        '[h,fig] = gcbo;'...
        'ud =  get(fig,''UserData'');'...
        'num = get(h,''UserData'');'...
        'ud.c.pname{num} = get(ud.h.pname(num),''string'');'...
        'ud.flag = 0;'...
        'set(gcf,''UserData'',ud);'];
    
    pvCallBack = [
        '[h,fig] = gcbo;'...
        'ud =  get(fig,''UserData'');'...
        'num = get(h,''UserData'');'...
        'ud.c.pval{num} = get(ud.h.pval(num),''string'');'...
        'ud.flag = 0;'...
        'set(gcf,''UserData'',ud);'];
    
    
    pnamew = 40*rr;
    %pvalw = 50*rr;
    peqw = 10*rr;
    pbot(3) = pfrbot + 5;
    pbot(2) = pbot(3) + separation;
    pbot(1) = pbot(2) + separation;
    
    %pleft1 = eqleft + 50*rr + 5;
    %peqleft1 = pleft1 + pnamew;
    %pvleft1 = peqleft1 + peqw;
    %pleft2 = dfigwidth - 10 - pnamew - pvalw - peqw;
    %peqleft2 = pleft2 + pnamew;
    %pvleft2 = peqleft2 + peqw;
    
    paratit=uicontrol('style','text',...
        'horizon','center',...
        'string',{'Parameters';'or';'expressions'},...
        'backgroundcolor',tcolor',...
        'visible','off');
    ext = get(paratit,'extent');
    paratitw = ext(3);
    pos = [eqleft pfrbot+2+texth/2 paratitw 2.1*texth];
    set(paratit,'pos',pos);
    psep = 20;
    pvalw = (dfigwidth - 2*eqleft - paratitw)/2 - psep - pnamew - peqw;
    pval = ud.c.pval;
    pname = ud.c.pname;
    for jj = 1:3
        for kk = 1:2
            pleft = eqleft + paratitw +psep +(kk-1)*(pnamew+peqw+pvalw+ psep);
            peqleft = pleft + pnamew;
            pvleft = peqleft + peqw;
            K = kk +2*(jj-1);
            name = pname{K};
            value = pval{K};
            ud.h.pname(K) = uicontrol('style','edit',...
                'pos',[pleft pbot(jj) pnamew texth],...
                'horizon','right','string',name,...
                'UserData',K,...
                'CallBack',pnCallBack,...
                'visible','off',...
                'backgroundcolor','w');
            [~] = uicontrol('style','text',...
                'pos',[peqleft pbot(jj)-fudge peqw texth],...
                'horizon','center',...
                'string','=',...
                'visible','off');
            
            ud.h.pval(K) = uicontrol('style','edit',...
                'pos',[pvleft pbot(jj) pvalw texth],...
                'string',value,...
                'CallBack',pvCallBack,...
                'visible','off',...
                'UserData',K,...
                'backgroundcolor','w');
        end
    end
    
    ud.c.pname = pname;
    ud.c.pval = pval;
    
    [~] = uicontrol('style','push',...
        'pos',qwind,...
        'string','Quit','CallBack',...
        'matpplane(''quit'')',...
        'visible','off');
    
    [~] = uicontrol('style','push',...
        'pos',rwind,...
        'string','Revert',...
        'CallBack','matpplane(''revert'')',...
        'visible','off');
    
    [~] = uicontrol('style','push',...
        'pos',pwind,...
        'string','Proceed',...
        'CallBack','matpplane(''proceed'')',...
        'visible','off');
    
    [~] = uicontrol('style','frame','pos',ffrwind,'visible','off');
    
    [~] = uicontrol('style','text',...
        'pos',[ffrleft+5,winbot1,ffrw-10,texth],...
        'string','The direction field.',...
        'horizon','center','visible','off');
    
    radleft = ffrleft + 3;
    radw = 50*rr;
    
    typewindw = radw +6;
    typewind = [ffrleft, ffrbot, typewindw, ffrht-separation-3];
    textwindl = ffrleft+typewindw;
    textleft = textwindl + 3;
    textw = ffrw - typewindw;
    textwind = [textwindl, ffrbot,textw, ffrht-separation-3];
    [~] = uicontrol('style','frame','pos',typewind,'visible','off');
    
    [~] = uicontrol('style','frame','pos',textwind,'visible','off');
    
    switch ud.o.fieldtype
        case 'nullclines'
            rval1 = 1;rval2 = 0;rval3 = 0;rval4 = 0;
        case 'lines'
            rval1 = 0;rval2 = 2;rval3 = 0;rval4 = 0;
        case 'arrows'
            rval1 = 0;rval2 = 0;rval3 = 3;rval4 = 0;
        case 'none'
            rval1 = 0;rval2 = 0;rval3 = 0;rval4 = 4;
        otherwise
            error('PPLANE:UnknownFieldType1',['Unknown fieldtype ',ud.o.fieldtype,'.'])
    end
    
    ud.h.rad(1) = uicontrol('style','radio',...
        'pos',[radleft winbot4 radw texth],...
        'string','Nullclines',...
        'value',rval1,...
        'visible','off');
    
    ud.h.rad(2) = uicontrol('style','radio',...
        'pos',[radleft winbot3 radw texth],...
        'string','Lines',...
        'value',rval2,...
        'max',2,...
        'visible','off');
    
    ud.h.rad(3) = uicontrol('style','radio',...
        'pos',[radleft winbot2 radw texth],...
        'string','Arrows',...
        'value',rval3,...
        'max',3,...
        'visible','off');
    
    ud.h.rad(4) = uicontrol('style','radio',...
        'pos',[radleft winbot5 radw texth],...
        'string','None',...
        'value',rval4,...
        'max',4,...
        'visible','off');
    
    for i=1:4
        set(ud.h.rad(i),'UserData',ud.h.rad(:,[1:(i-1),(i+1):4]));
    end
    
    CallBackrad = [
        'me = get(gcf,''currentobject'');',...
        'kk = get(me,''max'');',...
        'set(get(me,''UserData''),''value'',0),',...
        'set(me,''value'',kk);',...
        'ud = get(gcf,''UserData'');',...
        'switch kk,',...
        '       case 1, ud.c.fieldtype = ''nullclines'';',...
        '       case 2, ud.c.fieldtype = ''lines'';',...
        '       case 3, ud.c.fieldtype = ''arrows'';',...
        '       case 4, ud.c.fieldtype = ''none'';',...
        'end,',...
        'set(gcf,''UserData'',ud);'];
    
    set(ud.h.rad,'CallBack',CallBackrad);
    
    nfptsstr = {'Number of'; 'field points per'; 'row or column.'};
    [~] = uicontrol('style','text',...
        'pos',[textleft winbot4 textw-5 2.5*texth],...
        'string',nfptsstr,...
        'horizon','center',...
        'visible','off');
    
    CallBacknfpts = [
        'ppset = findobj(''name'',''matpplane Setup'');',...
        'ud = get(ppset,''UserData'');'...
        'me = ud.h.npts;',...
        'kk = str2num(get(me,''string''));',...
        'if isempty(kk),',...
        '  set(me,''string'',''?'');',...
        '  kk = NaN;',...
        'else,',...
        '  kk = floor(kk);',...
        '  [m,N] = computer;'...
        '  if (N <= 8192),',...
        '	N = 32;',...
        '  else,',...
        '	N = 50;',...
        '  end,'...
        '  kk = min([N,max([5,kk])]);'...
        '  set(me,''string'',num2str(kk));'...
        'end,'...
        'ud.c.npts = kk;',...
        'set(ppset,''UserData'',ud);'];
    
    npos = [textleft+(textw -30*rr)/2,winbot5 30*rr,texth];
    ud.h.npts = uicontrol('style','edit',...
        'pos',npos,...
        'string',ud.o.npts,...
        'CallBack',CallBacknfpts,...
        'backgroundcolor','w',...
        'visible','off');
    
    delgall = ['sud = get(gcf,''UserData'');',...
        'mh = get(sud.h.gallery,''children'');',...
        'add = findobj(sud.h.gallery,''tag'',''add system'');',...
        'mh(find(mh == add)) = [];',...
        'delete(mh);',...
        'set(sud.h.gallery,''UserData'',[]);',...
        'set(findobj(''tag'',''load default''),''enable'',''on'')'];
    
    
    % Menus
    
    hhsetup = get(0,'showhiddenhandles');
    set(0,'showhiddenhandles','on');
    mefile = findobj(ppset,'label','&File');
    meedit = findobj(ppset,'label','&Edit');
    delete(findobj(ppset,'label','&Tools'));
    delete(findobj(ppset,'label','&View'));
    delete(findobj(ppset,'label','&Insert'));
    
    % File menu
    
    meexp = findobj(mefile,'label','&Export...');
    meprev = findobj(mefile,'label','Print Pre&view...');
    mepset = findobj(mefile,'label','Pa&ge Setup...');
    set(get(mefile,'child'),'vis','off');
    [~] = uimenu(mefile,'label','Load a system ...',...
        'CallBack','matpplane(''loadsyst'',''system'');',...
        'pos',1);
    [~] = uimenu(mefile,'label','Save the current system ...',...
        'CallBack','matpplane(''savesyst'',''system'');',...
        'pos',2);
    [~] = uimenu(mefile,'label','Load a gallery ...',...
        'CallBack','matpplane(''loadsyst'',''gallery'');',...
        'separator','on','pos',3);
    [~] = uimenu(mefile,'label','Save a gallery ...',...
        'CallBack','matpplane(''savesyst'',''gallery'');',...
        'tag','savegal','pos',4);
    [~] = uimenu(mefile,'label','Delete the current gallery',...
        'CallBack',delgall,'pos',5);
    [~] = uimenu(mefile,'label','Load the default gallery',...
        'CallBack','matpplane(''loadsyst'',''default'');',...
        'enable','on',...
        'tag','load default','pos',6);
    [~] = uimenu(mefile,'label','Revert','CallBack',...
        'matpplane(''revert'')',...
        'separator','on','pos',7);
    [~] = uimenu(mefile,...
        'label','Proceed',...
        'CallBack','matpplane(''proceed'')',...
        'separator','off',...
        'accelerator','G','pos',8);
    set(mepset,'vis','on','pos',9);
    set(meprev,'vis','on','pos',10);
    set(meexp,'vis','on','pos',11,'separator','off');
    [~] = uimenu(mefile,'label',...
        'Restart matpplane',...
        'CallBack','matpplane(''restart'')',...
        'separator','on','pos',12);
    
    [~] = uimenu(mefile,...
        'label','Quit matpplane',...
        'CallBack','matpplane(''quit'')',...
        'separator','off','pos',13);
    
    % Edit menu
    
    set(get(meedit,'child'),'vis','off');
    [~] = uimenu(meedit,'label','Clear equations',...
        'CallBack',['ud = get(gcf,''UserData'');h = ud.h;',...
        'set([h.xvar,h.xder,h.yvar,h.yder],''string'','''');'],...
        'accelerator','E');
    
    pclear = [
        'ud = get(gcf,''UserData'');h = ud.h;',...
        'set([h.pname,h.pval],''string'','''');',...
        'ud.c.pname = {'''','''','''','''','''','''',''''};',...
        'ud.c.pval = {'''','''','''','''','''','''',''''};',...
        'set(gcf,''UserData'',ud);',...
        ];
    
    [~] = uimenu(meedit,'label','Clear parameters',...
        'CallBack',pclear,...
        'accelerator','N');
    
    [~] = uimenu(meedit,'label','Clear display window',...
        'CallBack',['ud = get(gcf,''UserData'');',...
        'set(ud.h.wind,''string'','''');'],...
        'accelerator','D');
    
    allclear = [
        'ud = get(gcf,''UserData'');h = ud.h;',...
        'set([h.xvar,h.xder,h.yvar,h.yder],''string'','''');',...
        'set([h.pname,h.pval,h.wind],''string'','''');',...
        'ud.c.pname = {'''','''','''','''','''','''',''''};',...
        'ud.c.pval = {'''','''','''','''','''','''',''''};',...
        'set(gcf,''UserData'',ud);',...
        ];
    
    [~] = uimenu(meedit,'label','Clear all',...
        'CallBack',allclear,...
        'accelerator','A',...
        'separator','on');
    
    % Gallery menu.
    
    sysmenu = uimenu('label','Gallery','visible','off','pos',3);
    
    [~] = uimenu(sysmenu,'label','Add current system to the gallery',...
        'CallBack','matpplane(''addgall'');','tag','add system');
    sep = 'on';
    for kk = 1:length(system)
        kkk = num2str(kk);
        if kk == 2, sep = 'off';end
        [~] = uimenu(sysmenu,'label',system(kk).name,...
            'CallBack',['matpplane(''system'',',kkk,')'],...
            'separator',sep,'visible','off');
    end
    set(sysmenu,'UserData',system);
    ud.h.gallery = sysmenu;
    ud.flag = 0;
    ud.egg = exist('EASTEREGG','file');
    
    % Record the handles in the User Data of the Set Up figure.
    
    set(ppset,'UserData',ud);
    hhhh = findobj(ppset,'type','uicontrol');
    set(hhhh,'units','normal')
    
    set(ppset,'visible','on','resize','on');
    set(get(ppset,'children'),'visible','on');
    set(get(sysmenu,'children'),'visible','on');
    %end
    set(0,'showhiddenhandles',hhsetup);
    
elseif strcmp(action,'savesyst')
    
    % Launch the Save System or Save Gallery dialog
    
    ppset = findobj('name','matpplane Setup');
    type = input1;
    sud = get(ppset,'UserData');
    
    switch type
        case 'system'
            %systems = get(sud.h.gallery,'UserData');
            
            newsyst = sud.c;
            fn = newsyst.name;
            if ~isempty(fn)
                %fn(find(abs(fn)==32))='_';   % Replace spaces by underlines.
                fn = strrep(fn,' ','_');
            end
            fn = [fn, '.pps'];
            %             comp = computer;
            %             switch  comp
            %                 case 'PCWIN'
            %                     filter = [sud.ppdir, '\',fn];
            %                 case 'MAC2'
            %                     filter = [sud.ppdir,':', fn];
            %                 otherwise
            %                     filter = [sud.ppdir,'/', fn];
            %             end
            filter = fullfile(sud.ppdir,fn);
            [fname,pname] = uiputfile(filter,'Save the system as:');
            if fname == 0,return;end
            if ~strcmp(fname,fn)
                ll = length(fname);
                if ll>4 && strcmp(fname(ll-3:ll),'.pps')
                    fn = fname;
                else
                    fn = [fname, '.pps'];
                end
                newsyst.name = fn;
                sud.c.name = fn;
                set(ppset,'UserData',sud);
            end
            newsysts = newsyst;
            
        case 'gallery'
            systems = get(sud.h.gallery,'UserData');
            ll = length(systems);
            if ll == 0
                warndlg('There are no systems to make up a gallery.','Warning');
                return
            end
            names = cell(ll,1);
            for j=1:ll
                names{j} = systems(j).name;
            end
            [sel,~] = listdlg('PromptString','Select the systems',...
                'Name','Gallery selection',...
                'ListString',names);
            if isempty(sel)
                return
            else
                newsysts = systems(sel);
            end
            %             comp = computer;
            %             switch  comp
            %                 case 'PCWIN'
            %                     prompt = [sud.ppdir,'\*.ppg'];
            %                 case 'MAC2'
            %                     prompt = [sud.ppdir,':*.ppg'];
            %                 otherwise
            %                     prompt = [sud.ppdir,'/*.ppg'];
            %             end
            prompt = fullfile(sud.ppdir,'*.ppg');
            [fname,pname] = uiputfile(prompt,'Save the gallery as:');
            ll = length(fname);
            if ll>4 && strcmp(fname(ll-3:ll),'.ppg')
                fn = fname;
            else
                fn = [fname, '.ppg'];
            end
            newsyst.name = fn(1:ll-4); %#ok<STRNU> % Why is this here?
            sud.c.name = fn(1:ll-4);
            set(ppset,'UserData',sud);
            
    end  % switch type
    
    ll = length(newsysts);
    fid = fopen([pname fn],'w');
    ppstring = '%%%% PPLANE file %%%%';
    fprintf(fid,[ppstring,'\n']);
    for k = 1:ll
        fprintf(fid,'\n');
        nstr = newsysts(k).name;
        nstr = strrep(nstr,'''','''''');
        nstr = ['H.name = ''', nstr, ''';\n']; %#ok<AGROW>
        fprintf(fid,nstr);
        xname = newsysts(k).xvar;
        xnstr = ['H.xvar = ''', xname];
        xnstr = strrep(xnstr,'\','\\');
        xnstr = [xnstr, ''';\n']; %#ok<AGROW>
        fprintf(fid,xnstr);
        yname = newsysts(k).yvar;
        ynstr = ['H.yvar = ''', yname];
        ynstr = strrep(ynstr,'\','\\');
        ynstr = [ynstr, ''';\n']; %#ok<AGROW>
        fprintf(fid,ynstr);
        xder = newsysts(k).xder;
        xdstr = ['H.xder = ''', xder];
        xdstr = strrep(xdstr,'\','\\');
        xdstr = [xdstr, ''';\n']; %#ok<AGROW>
        fprintf(fid,xdstr);
        yder = newsysts(k).yder;
        ydstr = ['H.yder = ''', yder];
        ydstr = strrep(ydstr,'\','\\');
        ydstr = [ydstr, ''';\n']; %#ok<AGROW>
        fprintf(fid,ydstr);
        
        pname = strrep(newsysts(k).pname,'\','\\');
        pval = strrep(newsysts(k).pval,'\','\\');
        pnl = length(pname);
        pvl = length(pval);
        for kk = 1:6
            if kk <= pnl
                pns = pname{kk};
            else
                pns = '';
            end
            if kk <= pvl
                pvs = pval{kk};
            else
                pvs = '';
            end
            if kk == 1
                pnstr = ['H.pname = {''', pns, ''''];
                pvstr = ['H.pval = {''', pvs, ''''];
            else
                pnstr = [pnstr, ',''',pns, '''']; %#ok<AGROW>
                pvstr = [pvstr, ',''',pvs, '''']; %#ok<AGROW>
            end
        end
        pnstr = [pnstr, '};\n']; %#ok<AGROW>
        pvstr = [pvstr, '};\n']; %#ok<AGROW>
        
        
        fprintf(fid,pnstr);
        fprintf(fid,pvstr);
        ftstr = ['H.fieldtype = ''', newsysts(k).fieldtype];
        ftstr = strrep(ftstr,'\','\\');
        ftstr = [ftstr, ''';\n']; %#ok<AGROW>
        fprintf(fid,ftstr);
        nstr = ['H.npts = ', num2str(newsysts(k).npts),';\n'];
        fprintf(fid,nstr);
        wind = newsysts(k).wind;
        wstr = ['H.wind = [', num2str(wind),'];\n'];
        fprintf(fid,wstr);
        
        
    end
    fclose(fid);
    
elseif strcmp(action,'loadsyst')
    
    %  This loads either a system or a gallery.
    
    sud = get(gcf,'UserData');
    pos = get(gcf,'pos');
    wpos = [pos(1),pos(2)+pos(4)+20,300,20];
    waith = figure('pos',wpos,...
        'vis','off',...
        'next','replace',...
        'menubar','none',...
        'resize','off',...
        'createfcn','');
    axes('pos',[0.01,0.01,0.98,0.98],...
        'vis','off');
    xp = [0 0 0 0];
    yp = [0 0 1 1];
    xl = [1 0 0 1 1];
    yl = [0 0 1 1 0];
    patchh = patch(xp,yp,'r','edgecolor','r','erase','none');
    lineh = line(xl,yl,'erase','none','color','k');
    type = input1;
    set(sud.h.gallery,'enable','off');
    if strcmp(type,'default')
        set(waith,'name','Loading the default gallery.','vis','on');
        set(findobj('tag','load default'),'enable','off');
        megall = sud.h.gallery;
        mh = get(megall,'children');
        add = findobj(megall,'tag','add system');
        mh(mh == add) = [];
        delete(mh);
        newsysstruct = get(megall,'UserData');
        system = sud.system;
        ll = length(system);
        x = 1/(ll+2);
        xp = [xp(2),x,x,xp(2)];
        set(patchh,'xdata',xp);
        set(lineh,'xdata',xl);
        drawnow;
        sep = 'on';
        for kk = 1:length(system)
            kkk = num2str(kk);
            if kk ==2, sep = 'off';end
            uimenu(megall,'label',system(kk).name,...
                'CallBack',['matpplane(''system'',',kkk,')'],...
                'separator',sep);
        end % for
        set(megall,'UserData',system);
    else
        %         comp = computer;
        %         switch  comp
        %             case 'PCWIN'
        %                 prompt = [sud.ppdir,'\'];
        %             case 'MAC2'
        %                 prompt = [sud.ppdir,':'];
        %             otherwise
        %                 prompt = [sud.ppdir,'/'];
        %         end
        
        if strcmp(type,'system')
            %prompt = [prompt,'*.pps'];
            prompt = fullfile(sud.ppdir,'*.pps');
            [fname,pname] = uigetfile(prompt,'Select a system to load.');
        elseif strcmp(type,'gallery')
            %prompt = [prompt,'*.ppg'];
            prompt = fullfile(sud.ppdir,'*.ppg');
            [fname,pname] = uigetfile(prompt,'Select a gallery to load.');
        end  % if strcmp
        
        if fname == 0
            delete(waith);
            set(sud.h.gallery,'enable','on');
            return;
        end
        set(waith,'name',['Loading ',fname],'vis','on');
        fid = fopen([pname fname],'r');
        sline = fgetl(fid);
        if strcmp(sline,'%% PPLANE file %%')
            the_date = 'new';
        else
            the_date = 'old';
        end
        newsysts = {};
        switch the_date
            case 'old'
                newsysts{1} = sline;
                kk = 1;
            case 'new'
                kk = 0;
        end
        while ~feof(fid)
            kk = kk + 1;
            newsysts{kk} = fgetl(fid); %#ok<AGROW>
        end
        fclose(fid);
        newsysts = newsysts(1:kk);
        falseFlag = 0;
        switch the_date
            case 'old'
                if mod(kk,19)
                    falseFlag = 1;
                end
            case 'new'
                if mod(kk,11)
                    falseFlag = 1;
                end
        end %switch date
        if falseFlag
            if strcmp(type,'system')
                warndlg(['The file ',fname, ' does not define a proper system.'],...
                    'Warning');
            elseif strcmp(type,'gallery')
                warndlg(['The file ',fname, ' does not define a proper gallery.'],...
                    'Warning');
            end
            set(sud.h.gallery,'enable','on');
            delete(waith);
            return
        end %if false
        switch the_date
            case 'old'
                x = 19/(kk+38);
                xp = [xp(2),x,x,xp(2)];
                set(patchh,'xdata',xp);
                set(lineh,'xdata',xl);
                drawnow;
                nnn = kk/19;
                %flds = fieldnames(sud.c);
                %flds=flds(:);
                for j = 0:(nnn-1)
                    newsystemp = newsysts((j*19+1):(j+1)*19);
                    newsyst.name = newsystemp{1};
                    newsyst.xvar = newsystemp{2};
                    newsyst.yvar = newsystemp{3};
                    newsyst.xder = newsystemp{4};
                    newsyst.yder = newsystemp{5};
                    newsyst.pname = {newsystemp{6}, newsystemp{7},...
                        newsystemp{8}, newsystemp{9},'',''};
                    newsyst.pval = {newsystemp{10},...
                        newsystemp{11},...
                        newsystemp{12},...
                        newsystemp{13},'',''};
                    newsyst.fieldtype = newsystemp{14};
                    newsyst.npts = str2double(newsystemp{15});
                    wind = newsystemp(16:19);
                    newsyst.wind = [str2double(wind{1}),str2double(wind{2}),...
                        str2double(wind{3}),str2double(wind{4})];
                    newsysstruct(j+1) = newsyst; %#ok<AGROW>
                end  % for j = 0:(nnn-1)
            case 'new'
                x = 11/(kk+22);
                xp = [xp(2),x,x,xp(2)];
                set(patchh,'xdata',xp);
                set(lineh,'xdata',xl);
                drawnow;
                nnn = kk/11;
                for j = 1:nnn
                    for k = 2:11
                        eval(newsysts{(j-1)*11+k});
                    end
                    newsysstruct(j) = H; %#ok<NODEF,AGROW>  % Somehow h is defined
                end
                
        end %switch date
    end  % if strcmp(type,'default') & else
    nnn = length(newsysstruct);
    ignoresyst = {};
    for j = 1:nnn
        x = (j+1)/(nnn+2);
        xp = [xp(2),x,x,xp(2)];
        set(patchh,'xdata',xp);
        set(lineh,'xdata',xl);
        drawnow;
        newsyst = newsysstruct(j);
        sname = newsyst.name;
        %sname(find(abs(sname) == 95)) = ' '; % Replace underscores with spaces.
        sname = strrep(sname,'_','');
        newsyst.name = sname;
        ignore = matpplane('addgall',newsyst);
        if ignore == -1
            ignoresyst{length(ignoresyst)+1} = sname; %#ok<AGROW>
        end
    end % for j = 1:nnn
    l = length(ignoresyst);
    if l  % There was at least one system which was a dup with a
        % different name.
        if l == 1
            message = {['The system ',ignoresyst{1},'" duplicates a ',...
                'system already in the gallery and was not added.']};
        else
            message = 'The systems ';
            for k = 1:(l-1)
                message = [message,'"',ignoresyst{k},'", ']; %#ok<AGROW>
            end
            message = {[message,'and "',ignoresyst{l},'" duplicate ',...
                'systems already in the gallery and were not added.']};
        end % if l == 1 & else
        helpdlg(message,'Ignored systems');
    end  % if l
    if strcmp(type,'system') % Added a system.
        if ignore > 0 % The system was ignored.
            kk = ignore;
        else
            systems = get(sud.h.gallery,'UserData');
            kk = length(systems);
        end
        matpplane('system',kk);
    end
    if strcmp('type','default')
        matpplane('system',1);
    end
    set(sud.h.gallery,'enable','on');
    x = 1;
    xp = [xp(2),x,x,xp(2)];
    set(patchh,'xdata',xp);
    set(lineh,'xdata',xl);
    drawnow;
    delete(waith);
    
elseif strcmp(action,'addgall')
    
    % Add all systems to the gallery
    
    output = 0;
    ppset = findobj('name','matpplane Setup');
    sud = get(ppset,'UserData');
    if nargin < 2    % We are adding the current system.
        
        syst = sud.c;
        snstr = 'Provide a name for this system.';
        sname = inputdlg(snstr,'System name',1,{syst.name});
        if isempty(sname),return;end
        sname = sname{1};
        if ~strcmp(sname,syst.name)
            sud.c.name = sname;
            set(ppset,'UserData',sud);
            syst.name = sname;
        end
        
    else  % We have a system coming from a file.
        syst = input1;
        sname = syst.name;
    end
    pnl = length(syst.pname);
    for kk = (pnl+1):6
        syst.pname{kk} = '';
    end
    pvl = length(syst.pval);
    for kk = (pvl+1):6
        syst.pval{kk} = '';
    end
    
    systems = get(sud.h.gallery,'UserData');
    ll = length(systems);
    kk = 1;
    while (kk<=ll) && (~strcmp(sname,systems(kk).name))
        kk = kk + 1;
    end
    nameflag = (kk<=ll);
    ssyst = rmfield(syst,'name');
    kk = 1;
    while (kk<=ll) && (~isequal(ssyst,rmfield(systems(kk),'name')))
        kk = kk + 1;
    end
    systflag = 2*(kk<=ll);
    flag = nameflag + systflag;
    switch flag
        case 1  % Same name but different system.
            mh = findobj(sud.h.gallery,'label',sname);
            prompt = {['The system "',sname,'", which you wish to ',...
                'add to the gallery has ',...
                'the same name as a different system ',...
                'already in the gallery.  Please ',...
                'specify the name for the newly added system.'],...
                'Specify the name for the old system.'};
            title = 'Two systems with the same name';
            lineno = 1;
            defans = {sname,sname};
            answer = inputdlg(prompt,title,lineno,defans);
            if isempty(answer),return,end
            sname = answer{1};
            systems(kk).name = answer{2};
            set(mh,'label',answer{2});
            output = kk;
        case 2  % Two names for the same system.
            oldname = systems(kk).name;
            mh = findobj(sud.h.gallery,'label',oldname);
            
            prompt = {['The system "',sname,'", which you wish to add ',...
                'to the gallery is the same as a system which is ',...
                'already in the gallery with the name "',oldname,'".  ',...
                'Please specify which name you wish to use.']};
            title = 'Two names for the same system.';
            lineno = 1;
            defans = {oldname};
            answer = inputdlg(prompt,title,lineno,defans);
            if isempty(answer),return,end
            systems(kk).name = answer{1};
            set(mh,'label',answer{1});
            output = kk;
        case 3 % Systems and names the same.
            output = -1;
        otherwise
    end  % switch
    set(sud.h.gallery,'UserData',systems);
    syst.name = sname;
    if flag <=1
        switch ll
            case 0
                systems = syst;
                sepstr = 'on';
            case 10
                systems(11) = syst;
                if strcmp(systems(10).name,'square limit set')
                    sepstr = 'on';
                else
                    sepstr = 'off';
                end
            otherwise
                systems(ll+1) = syst;
                sepstr = 'off';
        end
        kkk = num2str(ll+1);
        [~] = uimenu(sud.h.gallery,'label',sname,...
            'CallBack',['matpplane(''system'',',kkk,')'],...
            'separator',sepstr);
        set(findobj('tag','savegal'),'enable','on');
    end
    set(sud.h.gallery,'UserData',systems);
    
elseif strcmp(action,'system')
    
    % Launched when user clicks on a System in the Gallery menu of Setup Window
    
    ppset = findobj('name','matpplane Setup');
    ud = get(ppset,'UserData');
    kk = input1;
    if ischar(kk)
        kk = str2double(input1);
    end
    system = get(ud.h.gallery,'UserData');
    syst = system(kk);
    xname = syst.xvar;
    yname = syst.yvar;
    set(ud.h.xvar,'string',xname);
    set(ud.h.yvar,'string',yname);
    set(ud.h.xder,'string',syst.xder);
    set(ud.h.yder,'string',syst.yder);
    pname = syst.pname;
    pval = syst.pval;
    pnl = length(pname);
    pvl = length(pval);
    for kk = 1:6
        if kk <= pnl
            set(ud.h.pname(kk),'string',pname{kk});
        else
            set(ud.h.pname(kk),'string','');
            syst.pname{kk} = '';
        end
        if kk <= pvl
            set(ud.h.pval(kk),'string',pval{kk});
        else
            set(ud.h.pval(kk),'string','');
            syst.pval{kk} = '';
        end
    end
    ud.o = syst;
    ud.c = syst;
    set(ud.h.twind(1),'string',['The minimum value of ',xname,' = ']);
    set(ud.h.twind(2),'string',['The maximum value of ',xname,' = ']);
    set(ud.h.twind(3),'string',['The minimum value of ',yname,' = ']);
    set(ud.h.twind(4),'string',['The maximum value of ',yname,' = ']);
    for kk = 1:4
        set(ud.h.wind(kk),'string',num2str(syst.wind(kk)));
    end
    set(ud.h.npts,'string',num2str(syst.npts));
    
    switch syst.fieldtype
        case 'nullclines'
            rval = [1 0 0 0];
        case 'lines'
            rval = [0 2 0 0];
        case 'arrows'
            rval = [0 0 3 0];
        case 'none'
            rval = [0 0 0 4];
        otherwise
            error('PPLANE:UnknownFieldType2',['Unknown fieldtype ',ud.o.fieldtype,'.'])
    end
    
    for i=1:4
        set(ud.h.rad(i),'value',rval(i));
    end
    ud.flag = 0;
    set(ppset,'UserData',ud);
    
elseif strcmp(action,'revert')
    % When you click the "Revert" button, it is supposed to reset the data
    % in the setup window. Actually it does nothing
    % FIXME: This doesn't actually revert anything. ud.o and ud.c are the same
    ud = get(gcf,'UserData');
    ud.c = ud.o;
    syst = ud.o;
    xname = syst.xvar;
    yname = syst.yvar;
    set(ud.h.xvar,'string',xname);
    set(ud.h.yvar,'string',yname);
    set(ud.h.xder,'string',syst.xder);
    set(ud.h.yder,'string',syst.yder);
    pname = syst.pname;
    pval = syst.pval;
    pnl = length(pname);
    pvl = length(pval);
    for kk = 1:6
        if kk <= pnl
            nstr = pname(kk);
        else
            nstr = '';
        end
        if kk <= pvl
            vstr = pval(kk);
        else
            vstr = '';
        end
        set(ud.h.pname(kk),'string',nstr);
        set(ud.h.pval(kk),'string',vstr);
    end
    set(ud.h.twind(1),'string',['The minimum value of ',xname,' = ']);
    set(ud.h.twind(2),'string',['The maximum value of ',xname,' = ']);
    set(ud.h.twind(3),'string',['The minimum value of ',yname,' = ']);
    set(ud.h.twind(4),'string',['The maximum value of ',yname,' = ']);
    for kk = 1:4
        set(ud.h.wind(kk),'string',num2str(syst.wind(kk)));
    end
    set(ud.h.npts,'string',num2str(syst.npts));
    
    switch syst.fieldtype
        case 'lines'
            rval(1) = 1;
            rval(2) = 0;rval(3) = 0;
        case 'arrows'
            rval(1) = 0;rval(2) = 2;rval(3) = 0;
        case 'none'
            rval(1) = 0;rval(2) = 0;rval(3) = 3;
        otherwise
            error('PPLANE:UnknownFieldType1',['Unknown fieldtype ',ud.o.fieldtype,'.'])
    end
    
    for i=1:3
        set(ud.h.rad(i),'value',rval(i));
    end
    set(gcf,'UserData',ud);
    
elseif strcmp(action,'proceed')
    
    % Proceed connects Setup with the Display window.
    % Called when "Proceed" button pressed in setup window
    ppset = gcf;
    sud = get(ppset,'UserData');
    sud.o = sud.c;
    set(ppset,'UserData',sud);
    he = findobj('name','matpplane Equilibrium point data');
    hl = findobj('name','matpplane Linearization');
    close([he;hl]);
    % set([he;hl],'vis','off');
    
    % Some error checking that has to be done no matter what.
    
    WINvect = sud.c.wind;
    if any(isnan(WINvect))
        sud.flag = 0;
        set(ppset,'UserData',sud);
        errmsg = ['One of the entries defining the display window ',...
            'is not a number.'];
        fprintf('\a')
        errordlg(errmsg,'MATPPLANE error','on');
        return
    end
    xstr = sud.c.xvar;
    if isempty(xstr)
        sud.flag = 0;
        set(ppset,'UserData',sud);
        errmsg = 'The first dependent variable needs a name.';
        fprintf('\a')
        errordlg(errmsg,'MATPPLANE error','on');
        return
    end
    ystr = sud.c.yvar;
    if isempty(ystr)
        sud.flag = 0;
        set(ppset,'UserData',sud);
        errmsg = 'The second dependent variable needs a name.';
        fprintf('\a')
        errordlg(errmsg,'MATPPLANE error','on');
        return
    end
    if WINvect(2)<= WINvect(1)
        sud.flag = 0;
        set(ppset,'UserData',sud);
        errmsg = ['The minimum value of ', xstr,...
            ' must be smaller than the maximum value.'];
        fprintf('\a')
        errordlg(errmsg,'MATPPLANE error','on');
        return
    end
    if WINvect(4)<= WINvect(3)
        sud.flag = 0;
        set(ppset,'UserData',sud);
        errmsg = ['The minimum value of ', ystr,...
            ' must be smaller than the maximum value.'];
        fprintf('\a')
        errordlg(errmsg,'MATPPLANE error','on');
        return
    end
    if isnan(sud.c.npts)
        sud.flag = 0;
        set(ppset,'UserData',sud);
        errmsg = 'The entry for the number of field points is not a number.';
        fprintf('\a')
        errordlg(errmsg,'MATPPLANE error','on');
        return
    end
    
    %  sud.flag = 0 if this is the first time through for this equation,
    %  sud.flag = 1 if only the window dimensions or the field data
    %  have been changed.
    
    % If sud.flag == 1 we only have to update things.
    
    if (sud.flag == 1)
        Arrflag = sud.c.fieldtype;
        NumbFPts = sud.c.npts;
        ppdisp = findobj('name','matpplane Display');
        dud = get(ppdisp,'UserData');
        aud = get(dud.axes,'UserData');
        wind = sud.c.wind(:);
        if (~all(wind == dud.syst.wind(:)))
            dwind = [wind(1); wind(3); -wind(2); -wind(4)];
            aud.DY = [wind(2)-wind(1); wind(4)-wind(3)];
            aud.cwind = dwind - dud.settings.magn*[aud.DY;aud.DY];
            set(dud.axes,'UserData',aud);
        end
        
        arr = dud.arr;
        menull = findobj('tag','null');
        switch Arrflag
            case 'nullclines'
                %       set([arr.hx;arr.hy;arr.barrows],'vis','on');
                set([arr.hx;arr.hy],'vis','on');
                set([arr.lines;arr.arrows],'vis','off');
                set(menull,'enable','on','label','Hide nullclines.');
            case 'lines'
                set(arr.lines,'vis','on');
                set([arr.hx;arr.hy;arr.arrows;arr.barrows],'vis','off');
                set(menull,'enable','on','label','Show nullclines.');
            case 'arrows'
                set(arr.arrows,'vis','on');
                set([arr.lines;arr.hx;arr.hy;arr.barrows],'vis','off');
                set(menull,'enable','on','label','Show nullclines.');
            otherwise
                set([arr.hx;arr.hy;arr.lines;arr.arrows;arr.barrows],'vis','off');
                set(menull,'enable','on','label','Show nullclines.');
        end
        dud.syst.fieldtype = Arrflag;
        set(ppdisp,'UserData',dud);
        
        if ( (NumbFPts ~= dud.syst.npts) || (any(WINvect ~= dud.syst.wind) ) )
            dud.syst.wind = WINvect;
            dud.syst.npts = NumbFPts;
            set(ppdisp,'UserData',dud);
            matpplane('dirfield',ppdisp);
        end
        figure(ppdisp);
    else
        sud.flag = 1;
        set(ppset,'UserData',sud);
        
        sud = get(ppset,'UserData');
        %  WINvect = sud.c.wind;
        
        Xname = sud.c.xvar;
        Yname = sud.c.yvar;
        xderivstr = sud.c.xder;
        yderivstr = sud.c.yder;
        pname = sud.c.pname;
        parav = sud.c.pval;
        
        % Convert the parameters to their current values.  First remove the
        % blanks.  Also remove the periods inserted by UserDatas attempting to
        % make the function array smart.
        
        %xderivstr(find(abs(xderivstr)==32))=[];
        xderivstr = strrep(xderivstr,' ','');
        l=length(xderivstr);
        for  k = fliplr(strfind(xderivstr,'.'))
            if (find('*/^' == xderivstr(k+1)))
                xderivstr = [xderivstr(1:k-1), xderivstr(k+1:l)];
            end
            l=l-1;
        end
        
        %yderivstr(find(abs(yderivstr)==32))=[];
        yderivstr = strrep(yderivstr,' ','');
        l=length(yderivstr);
        for k = fliplr(strfind(yderivstr,'.'))
            if (find('*/^' == yderivstr(k+1)))
                yderivstr = [yderivstr(1:k-1), yderivstr(k+1:l)];
            end
            l=l-1;
        end
        
        for kk = 1:6
            pval = parav{kk};
            if ~isempty(pval)
                %pabs = abs(pval);
                %kkk = find(pabs==32);
                %pval(kkk) = [];
                pval = strrep(pval,' ','');
                l = length(pval);
                for k = fliplr(strfind(pval,'.'))
                    if (find('*/^' == pval(k+1)))
                        pval = [pval(1:k-1), pval(k+1:l)];
                    end
                    l=l-1;
                end
                parav{kk} = pval;
            end
        end
        
        % Build strings for the title.
        txderstr = xderivstr;
        tyderstr = yderivstr;
        
        kxder = find(abs(txderstr)==42);
        txderstr(kxder)=' '*ones(size(kxder));
        txderstr = strrep(txderstr,'-',' - ');
        txderstr = strrep(txderstr,'+',' + ');
        %        if (abs(txderstr(1)) == 32)
        if strcmp(txderstr(1),' ')
            txderstr = txderstr(2:length(txderstr));
        end
        
        kxder = find(abs(tyderstr)==42);
        tyderstr(kxder)=' '*ones(size(kxder));
        tyderstr = strrep(tyderstr,'-',' - ');
        tyderstr = strrep(tyderstr,'+',' + ');
        if strcmp(tyderstr(1),' ')
            %if (abs(tyderstr(1)) == 32)
            tyderstr = tyderstr(2:length(tyderstr));
        end
        
        tstr1 = [Xname,' '' = ', txderstr];
        tstr2 = [Yname,' '' = ', tyderstr];
        tstr = char(tstr1,tstr2);
        dud.tstr = tstr; %#ok<STRNU>  % This seems to be needed
        pstr1 = {' ';' '};
        pstr2 = {' ';' '};
        pstr3 = {' ';' '};
        
        pstring = cell(6,1);
        for kk = 1:6
            if ~isempty(parav{kk})
                tpstr = parav{kk};
                %kxder = find(abs(tpstr)==42);  % Get rid of *s
                kxder = strfind(tpstr,'*');
                tpstr(kxder)=' '*ones(size(kxder));
                tpstr = strrep(tpstr,'-',' - ');  % Extra spaces
                tpstr = strrep(tpstr,'+',' + ');
                %if (abs(tpstr(1)) == 32)  % Get rid of starting space
                if strcmp(tpstr(1),' ')
                    tpstr = tpstr(2:length(tpstr));
                end
                pstring{kk} = [pname{kk},' = ', tpstr];
            else
                % pstring{kk} = [pname{kk},' = ',parav{kk}];
                pstring{kk} = [pname{kk},' = '];
            end
        end
        
        % Get ready to do some error trapping.
        
        SS = warning;
        warning off
        XxXxXx = WINvect(1) + rand*(WINvect(2)-WINvect(1));
        YyYyYy = WINvect(3) + rand*(WINvect(4)-WINvect(3));
%        err = 0;
        
        % Now we remove the backslashes (\) used to get Greek into the
        % labels.
        txname = Xname;
        tyname = Yname;
        xderivstr = strrep(xderivstr,'\',[]);
        yderivstr = strrep(yderivstr,'\',[]);
        Xname = strrep(Xname,'\',[]);
        Yname = strrep(Yname,'\',[]);

        
%        eval([Xname,'=XxXxXx;'],'err = 1;');
%         if err
%             sud.flag = 0;
%             set(ppset,'UserData',sud);
%             errmsg = ['"',xstr, '" is not a valid variable name in MATLAB.'];
%             fprintf('\a')
%             errordlg(errmsg,'MATPPLANE error','on');
%             return
%         end
%         err = 0;
%        eval([Yname,'=YyYyYy;'],'err = 2;');
%         if err
%             sud.flag = 0;
%             set(ppset,'UserData',sud);
%             errmsg = ['"',ystr, '" is not a valid variable name in MATLAB.'];
%             fprintf('\a')
%             errordlg(errmsg,'MATPPLANE error','on');
%             return
%         end
        
        if ~isvarname(Xname)
            sud.flag = 0;
            set(ppset,'UserData',sud);
            errmsg = ['"',xstr, '" is not a valid variable name in MATLAB.'];
            fprintf('\a')
            errordlg(errmsg,'MATPPLANE error','on');
            return
        end
        
        
        if ~isvarname(Yname)
            sud.flag = 0;
            set(ppset,'UserData',sud);
            errmsg = ['"',ystr, '" is not a valid variable name in MATLAB.'];
            fprintf('\a')
            errordlg(errmsg,'MATPPLANE error','on');
            return
        end        
        
        % Replace the parameters/expressions with their values.
        
        pflag = zeros(1,6);
        perr = [];
        for kk = 1:6
            if ~isempty(pname{kk})
                pn = pname{kk};
                pv = parav{kk};
                if isempty(pv)
                    perr = [perr, sud.h.pval(kk)]; %#ok<AGROW>
                else
                    if isempty(str2double(pv)) % This is an expression.
                        %tpv = pv;
                        %tpv(find(abs(tpv)==92))=[];
                        tpv=strrep(pv,'\','');
                        err = 0; pval = '';
                        eval(['pval = ',tpv,';'],'err=1;');
                        if err || isempty(pval)
                            errmsg = ['The expression for ',pn,' is not valid.'];
                            fprintf('\a')
                            errordlg(errmsg,'matpplane error','on');
                            warning(SS)
                            return
                        end
                    end
                    xxderivstr = matpplane('paraeval',pn,pv,xderivstr);
                    yyderivstr = matpplane('paraeval',pn,pv,yderivstr);
                    if ~strcmp(xxderivstr,xderivstr) || ~strcmp(yyderivstr,yderivstr)
                        pflag(kk) = 1;
                        xderivstr = xxderivstr;
                        yderivstr = yyderivstr;
                    end
                end
            end
        end
        
        % We have to make the derivative strings array smart.
        
        l = length(xderivstr);
        for k=fliplr(find((xderivstr=='^')|(xderivstr=='*')|(xderivstr=='/')))
            xderivstr = [xderivstr(1:k-1) '.' xderivstr(k:l)];
            l = l+1;
        end
        
        l = length(yderivstr);
        for k=fliplr(find((yderivstr=='^')|(yderivstr=='*')|(yderivstr=='/')))
            yderivstr = [yderivstr(1:k-1) '.' yderivstr(k:l)];
            l = l+1;
        end
        
        % Some more error trapping.
        
%        err = 0;res = 1;
%        eval(['res = ',xderivstr, ';'],'err = 1;');
%        if err || isempty(res)
         xderivfcn=str2func(['@(' Xname ',' Yname ')' xderivstr]);
         try
             xderivfcn(XxXxXx,YyYyYy);
         catch 
            if isempty(perr)
                errmsg = ['The first differential equation ',...
                    'is not entered correctly.'];
            else
                errstr1 = ['The first differential equation ',...
                    'does not evaluate correctly.'];
                errstr2 = ['At least one of the parameter values is not ',...
                    'a number.'];
                errmsg = char(errstr1,errstr2);
                
                set(perr,'string','?');
            end
            sud.flag = 0;
            set(ppset,'UserData',sud);
            fprintf('\a')
            errordlg(errmsg,'MATPPLANE error','on');
            return;
        end
        
%         err = 0;res = 1;
%         eval(['res = ',yderivstr, ';'],'err = 1;');
%         if err || isempty(res)
         yderivfcn=str2func(['@(' Xname ',' Yname ')' yderivstr]);
         try
             yderivfcn(XxXxXx,YyYyYy);
         catch 

            if isempty(perr)
                errmsg = ['The second differential equation ',...
                    'is not entered correctly.'];
            else
                errstr1 = ['The second differential equation ',...
                    'does not evaluate correctly.'];
                errstr2 = ['At least one of the parameter values is not ',...
                    'a number.'];
                errmsg = char(errstr1,errstr2);
                set(perr,'string','?');
            end
            sud.flag = 0;
            set(ppset,'UserData',sud);
            fprintf('\a')
            errordlg(errmsg,'MATPPLANE error','on');
            return;
        end
        
        % If an old function m-file exists delete it, and then build a new one.
        
        %    if (~strcmp(dfcn,'') & exist(dfcn)==2) delete([dfcn,'.m']);end
        tee = clock;
        tee = ceil(tee(6)*100);
        dfcn=['pptp',num2str(tee)];
        fcnstr = ['function YYyYypr = ',dfcn,'(t,YyYy)\n\n'];
        commstr = '%%%% Created by matpplane\n\n';
        varstr = [Xname,' = YyYy(1,:);', Yname,' = YyYy(2,:);\n\n'];
        lenstr = 'l = length(YyYy(1,:));\n';
        derstr1 = ['XxXxxpr = ', xderivstr,';\n'];
        derstr2 = 'if (length(XxXxxpr) < l) XxXxxpr = XxXxxpr*ones(1,l);end\n';
        derstr3 = ['YyYyypr = ', yderivstr,';\n'];
        derstr4 = 'if (length(YyYyypr) < l) YyYyypr = YyYyypr*ones(1,l);end\n';
        derstr5 = 'YYyYypr = [XxXxxpr;YyYyypr];\n';
        ppf = fopen([tempdir,dfcn,'.m'],'w');
        
        fprintf(ppf,fcnstr);
        fprintf(ppf,commstr);
        fprintf(ppf,varstr);
        fprintf(ppf,lenstr);
        fprintf(ppf,derstr1);
        fprintf(ppf,derstr2);
        fprintf(ppf,derstr3);
        fprintf(ppf,derstr4);
        fprintf(ppf,derstr5);
        fclose(ppf);
        
        
        
        % Find matpplane Display if it exists.
        % If matpplane Display exists, update it.  If it does not build it.
        
        ppdisp = findobj('name','matpplane Display');
        if (~isempty(ppdisp))
            figure(ppdisp);
            dud = get(ppdisp,'UserData');
            dud.syst = sud.c;
            dud.settings = sud.settings;
            dfcnn = dud.function;
            if  exist(dfcnn,'file')
                delete([tempdir,dfcnn,'.m']);
            end
            xmstr = [txname,' vs. t'];
            ymstr = [tyname,' vs. t'];
            set(dud.menu(3),'label',xmstr);
            set(dud.menu(5),'label',ymstr);
            menull = findobj('tag','null');
            if ~isempty(menull)
                delete(get(menull,'UserData'));
                set(menull,'UserData',[]);
            end
        else
            ppdisp = figure('name','matpplane Display',...
                'interrupt','on',...
                'visible','off',...
                'tag','matpplane');
            matpplane('figdefault',ppdisp);
            dud = get(ppdisp,'UserData');
            dud.syst = sud.c;
            switch dud.syst.name
                case 'pendulum'
                    dud.level = 'omega^2 - 2*cos(theta)';
                otherwise
                    dud.level = ' ';
            end
            dud.settings = sud.settings;
            dud.egg = sud.egg;
            dud.noticeflag = 1;
            dud.contours = zeros(0,1);
            fs = dud.fontsize;
            ssize = dud.ssize;
            r = ssize/10;
            ppaxw = 437*1.2;	% Default axes width
            ppaxh = 315*1.2;	% Default axes height
            ppaxl = 45*1.2;	% Default axes left
            buttw = 40*1.2;     % Default button width
            titleh = 45;    % Default title height.  This is changed later.
            nframeh = 70;   % Default notice frame height
            %ppaxb = 4+nframeh+35;
            bottomedge = 38;
            ppdh = bottomedge + nframeh + ppaxh + titleh;
            uni = get(0,'units');
            set(0,'units','pixels');
            ss = get(0,'screensize');
            set(0,'units',uni);
            sw = ss(3);sh = ss(4);
            bottom = 10;
            if r*ppdh > sh -bottom -35
                r = (sh-bottom-35)/ppdh;
                fs = 10*r;
                lw = 0.5*r;
                set(gcf,'defaultaxesfontsize',fs,'defaultaxeslinewidth',lw);
                set(gcf,'defaulttextfontsize',fs);
                set(gcf,'defaultlinelinewidth',lw);
                set(gcf,'defaultuicontrolfontsize',fs*0.9);
            end
            
            % Set up the bulletin window.
            
            nframe = uicontrol('style','frame','visible','on');
            nstr = {'More';'than';'five';'lines';'of text'};
            dud.notice = uicontrol('style','text',...
                'horiz','left',...
                'string',nstr,'visible','on');
            ext = get(dud.notice,'extent');
            nframeh = ext(4)+2;
            
            %titleh = r*titleh;
            ppaxl = r*ppaxl;
            ppaxw = r*ppaxw;
            ppaxh = r*ppaxh;
            ppaxb = nframeh+r*bottomedge;
            buttw = r*buttw;
            butth = fs+10*r;
            buttl = ppaxl + ppaxw + 5;
            buttsep = (ppaxh - butth)/2;
            % Set up the coordinate display
            
            cstr = '(0.99999999, 0,99999999)';
            dud.ccwind = uicontrol('style','text',...
                'horiz','left',...
                'string',cstr,...
                'visible','on');
            cext = get(dud.ccwind,'extent');
            ccwindtxt = uicontrol('style','text',...
                'horiz','left',...
                'string','Cursor position: ',...
                'visible','on');
            cwh = cext(4);
            cww = cext(3);
            
            
            % Set up the plot axes.
            
            dud.axes = axes('units','pix',...
                'position',[ppaxl,ppaxb,ppaxw,ppaxh],...
                'next','add',...
                'box','on',...
                'interrupt','on',...
                'xgrid','on',...
                'ygrid','on',...
                'drawmode','fast',...
                'visible','off',...
                'tag','display axes');
            
            % Set up the title.
            
            dud.title.axes = axes('box','off','xlim',[0 1],'ylim',[0 1],...
                'units','pix','vis','off',...
                'xtick',-1,'ytick',-1,...
                'xticklabel','','yticklabel','');
            
            dud.title.eq = text(0.01,0.5,' ','vert','middle');
            dud.title.p1 = text(0.75,0.5,' ','vert','middle');
            dud.title.p2 = text(0.65,0.5,' ','vert','middle');
            dud.title.p3 = text(0.55,0.5,' ','vert','middle');
            tstr = {'x_2';'y^2'};
            set(dud.title.eq,'string',tstr,'units','pix');
            ext = get(dud.title.eq ,'extent');
            titleh = ext(4)+15*r;
            set(dud.title.eq,'units','data');
            taxpos = [ppaxl,ppaxb+ppaxh,ppaxw,titleh];
            set(dud.title.axes,'pos',taxpos,'color',get(gcf,'color'));
            
            % Finish the positions.
            
            ppdw = buttl + buttw +5;
            ppdh = ppaxb+ppaxh+titleh;
            set(nframe,'pos',[10,1,ppdw-20,nframeh]);
            set(dud.notice,'pos',[15,1,ppdw-30,nframeh-2],...
                'string',{' ';' ';' ';' ';' '});
            
            ctext = get(ccwindtxt,'extent');
            cc1pos = [ppaxl,2+nframeh,ctext(3),cwh];
            cc2pos = [ppaxl+ctext(3),2+nframeh,cww,cwh];
            set(ccwindtxt,'pos',cc1pos);
            set(dud.ccwind,'pos',cc2pos,...
                'string','');
            
            ppdleft = max((sw-ppdw)/2,sw-ppdw-60);
            ppdbot = sh - ppdh - 35;
            ppdpos = [ppdleft,ppdbot,ppdw,ppdh];
            set(ppdisp,'resize','on');
            set(ppdisp,'pos',ppdpos);
            
            %Arrflag = sud.c.fieldtype;
            
            % Set up the buttons
            
            stopstr = 'aud = get(gca,''UserData'');aud.stop = 4;set(gca,''UserData'',aud);';
            
            dbutt(1) = uicontrol('style','push',...
                'pos',[buttl,ppaxb+2*buttsep,buttw,butth],...
                'string','Stop','CallBack',stopstr,...
                'vis','off','tag','stop');
            
            dbutt(2) = uicontrol('style','push',...
                'pos',[buttl,ppaxb,buttw,butth],...
                'string','Quit',...
                'CallBack','matpplane(''quit'')','visible','off');
            
            dbutt(3) = uicontrol('style','push',...
                'pos',[buttl,ppaxb+buttsep,buttw,butth],...
                'string','Print',...
                'CallBack','matpplane(''print'')','visible','off');
            
            dud.butt = dbutt;
            
            % Menus and Toolbar
            
            hhsetup = get(0,'showhiddenhandles');
            set(0,'showhiddenhandles','on');
            
            % Configure the Toolbar.
            
            %fixtb = 'set(gcbo,''state'',''off'');';
            
            set(ppdisp,'ToolBar','none');
            
            % Menus
            
            tmenu = findobj(ppdisp,'label','&Tools');
            delete(tmenu);
            
            % Insert menu
            
            imenu = findobj(gcf,'label','&Insert');
            inschild = get(imenu,'child');
            legitem = findobj(inschild,'label','&Legend');
            colitem = findobj(inschild,'label','&Colorbar');
            delete([legitem,colitem]);
            
            % File menu
            
            fmenu = findobj(ppdisp,'label','&File',...
                'pos',1);
            delete(findobj(fmenu,'label','&New Figure'));
            delete(findobj(fmenu,'label','&Open...'));
            delete(findobj(fmenu,'label','&Close'));
            set(findobj(fmenu,'label','&Save'),...
                'pos',1,'separator','off');
            set(findobj(fmenu,'label','Save &As...'),...
                'pos',2);
            set(findobj(fmenu,'label','&Export...'),...
                'pos',3);
            delete(findobj(fmenu,'label','Pre&ferences...'));
            set(findobj(fmenu,'label','Pa&ge Setup...'),'pos',4);
            set(findobj(fmenu,'label','Print Set&up...'),'pos',5);
            set(findobj(fmenu,'label','Print Pre&view...'),'pos',6);
            set(findobj(fmenu,'label','&Print...'),'pos',7);
            [~] = uimenu(fmenu,'label',...
                'Restart matpplane',...
                'CallBack','matpplane(''restart'')',...
                'separator','on');
            [~] = uimenu(fmenu,'label','Quit matpplane',...
                'CallBack','matpplane(''quit'')','separator','off');
            
            % Edit menu
            
            emenu = findobj(ppdisp,'label','&Edit',...
                'pos',2);
            menu(2) = uimenu(emenu,'label','Zoom in.',...
                'CallBack','matpplane(''zoomin'')',...
                'pos',1);
            
            [~] = uimenu(emenu,'label','Zoom in square.',...
                'CallBack','matpplane(''zoominsq'')',...
                'pos',2);
            
            [~] = uimenu(emenu,'label','Zoom back.','CallBack',...
                'matpplane(''zoomback'')',...
                'enable','on','tag','zbmenu',...
                'pos',3);
            
            [~] = uimenu(emenu,'label','Erase all solutions.',...
                'CallBack','matpplane(''dallsol'')',...
                'separator','on','pos',4);
            
            [~] = uimenu(emenu,'label','Erase all equilibrium points.',...
                'CallBack','matpplane(''dallep'')',...
                'separator','off',...
                'pos',5);
            
            [~] = uimenu(emenu,'label','Erase all marked initial points.',...
                'CallBack','matpplane(''dallics'')',...
                'separator','off',...
                'pos',6);
            
            [~] = uimenu(emenu,'label','Erase all level curves.',...
                'CallBack','matpplane(''dalllev'')',...
                'separator','off',...
                'pos',7);
            
            [~] = uimenu(emenu,'label','Erase all graphics objects.',...
                'CallBack','matpplane(''dall'')',...
                'separator','off',...
                'pos',8);
            
            [~] = uimenu(emenu,'label','Delete a graphics object.',...
                'CallBack','matpplane(''delete'')',...
                'visible','on',...
                'pos',9);
            
            [~] = uimenu(emenu,...
                'label','Enter text on the Display Window.',...
                'CallBack','matpplane(''text'')',...
                'separator','on',...
                'pos',10);
            set(findobj(emenu,'label','&Undo'),'separator','on',...
                'pos',11);
            set(findobj(emenu,'label','Cu&t'),'pos',12);
            set(findobj(emenu,'label','&Copy'),'pos',13);
            set(findobj(emenu,'label','&Paste'),'pos',14);
            set(findobj(emenu,'label','Clea&r'),'pos',15);
            set(findobj(emenu,'label','&Select All'),'pos',16);
            set(findobj(emenu,'label','Copy &Figure'),'pos',17);
            set(findobj(emenu,'label','Copy &Options'),'pos',18);
            set(findobj(emenu,'label','F&igure Properties'),'pos',19);
            set(findobj(emenu,'label','&Axes Properties'),'pos',20);
            set(findobj(emenu,'label','C&urrent Object Properties'),'pos',21);
            
            % Graph menu
            
            megraph  = uimenu('label','Graph','visible','off','pos',4);
            menu(3) = uimenu(megraph,'label',[txname,' vs. t'],...
                'CallBack','matpplane(''plotxy'',1)');
            
            menu(5) = uimenu(megraph,'label',[tyname,' vs. t'],...
                'CallBack','matpplane(''plotxy'',2)');
            
            [~] = uimenu(megraph,'label','Both',...
                'CallBack','matpplane(''plotxy'',3)');
            
            [~] = uimenu(megraph,'label','3 D',...
                'CallBack','matpplane(''plotxy'',4)');
            
            [~] = uimenu(megraph,'label','Composite',...
                'CallBack','matpplane(''plotxy'',5)');
            
            % Solutions menu
            
            solmenu = uimenu('label','Solutions','pos',3);
            [~] = uimenu(solmenu,'label','Keyboard input.','CallBack',...
                'matpplane(''kbd'')');
            
            [~]   = uimenu(solmenu,'label','Plot several solutions.',...
                'CallBack','matpplane(''several'')');
            
            [~]  = uimenu(solmenu,'label','Find an equilibrium point.',...
                'CallBack','matpplane(''eqpt'')','separator','on');
            
            menu(4) = uimenu(solmenu,...
                'label','List computed equilibrium points.',...
                'CallBack','matpplane(''eqptlist'')');
            
            [~]= uimenu(solmenu,...
                'label','Plot stable and unstable orbits.',...
                'CallBack','matpplane(''stunst'')','interrupt','on');
            
            meperiod = uimenu(solmenu,...
                'label', 'Find a nearly closed orbit');
            periodstr = ['ud = get(gcf,''UserData'');',...
                'me = gcbo;',...
                'ud.dir = get(me,''UserData'');',...
                'set(gcf,''UserData'',ud);',...
                'matpplane(''periodic'');'];
            dud.period(1) = uimenu(meperiod,...
                'label','forward',...
                'UserData',1,...
                'CallBack', periodstr);
            dud.period(2) = uimenu(meperiod,...
                'label','backward',...
                'UserData',-1,...
                'CallBack', periodstr);
            dud.period(3) = uimenu(meperiod,...
                'label','in both directions',...
                'UserData',0,...
                'CallBack', periodstr);
            
            nullCallBack = ['me = gcbo;',...
                'dud = get(gcf,''UserData'');',...
                'handx = [dud.arr.hx; dud.arr.hy];',...
                'handr = [dud.arr.hr; dud.arr.hth];',...
                'arron = get(dud.arr.arrows,''vis'');',...
                'lab = get(me,''label'');',...
                'switch lab,',...
                '  case ''Show nullclines.'',',...
                '    set(handx,''vis'',''on'');',...
                '    if strcmp(arron,''on''),',...
                '      set(dud.arr.barrows,''vis'',''off'');',...
                '    else,',...
                '      set(dud.arr.barrows,''vis'',''on'');',...
                '    end,',...
                '    set(handr,''vis'',''off'');',...
                '    set(me,''label'',''Hide nullclines.'');',...
                '    rnh = findobj(''tag'',''rnull'');'...
                '    set(rnh,''label'',''Show polar nullclines.'');',...
                '  case ''Hide nullclines.'','...
                '    set(handx,''vis'',''off'');',...
                '    set(dud.arr.barrows,''vis'',''off'');',...
                '    set(me,''label'',''Show nullclines.'');',...
                'end'];
            
            rnullCallBack = ['me = gcbo;',...
                'dud = get(gcf,''UserData'');',...
                'handx = [dud.arr.hx; dud.arr.hy];',...
                'handr = [dud.arr.hr; dud.arr.hth];',...
                'lab = get(me,''label'');',...
                'switch lab,',...
                '  case ''Show polar nullclines.'',',...
                '    set(handr,''vis'',''on'');',...
                '    set(handx,''vis'',''off'');',...
                '    set(dud.arr.barrows,''vis'',''off'');',...
                '    set(me,''label'',''Hide polar nullclines.'');',...
                '    nh = findobj(''tag'',''null'');'...
                '    set(nh,''label'',''Show nullclines.'');',...
                '  case ''Hide polar nullclines.'','...
                '    set(handr,''vis'',''off'');',...
                '    set(me,''label'',''Show polar nullclines.'');',...
                'end'];
            
            [~] = uimenu(solmenu,'label','Show nullclines.',...
                'CallBack',nullCallBack,'separator','on','tag','null');
            
            [~] = uimenu(solmenu,'label','Plot level curves.',...
                'CallBack','matpplane(''level'')',...
                'separator','off','tag','level');
            
            if dud.egg
                [~] = uimenu(solmenu,'label','Show polar nullclines.',...
                    'CallBack',rnullCallBack,'separator','off','tag','rnull');
                [~] = uimenu(solmenu,'label','Test case',...
                    'CallBack','matpplane(''test case'')',...
                    'separator','on','tag','testcase');
                
            end
            
            % Options menu
            
            menu(1) = uimenu('label','Options','visible','off','pos',5);
            [~]   = uimenu(menu(1),'label','Settings.',...
                'CallBack','matpplane(''settings'')');
            
            mesolve = uimenu(menu(1),'label','Solver.');
            
            solverstr = ['ud = get(gcf,''UserData'');',...
                'me = gcbo;',...
                'meud = get(me,''UserData'');',...
                'ud.settings.refine = meud.refine;',...
                'ud.settings.tol = meud.tol;',...
                'ud.settings.solver = meud.solver;',...
                'ud.settings.stepsize = meud.stepsize;',...
                'set(ud.solver,''checked'',''off'');',...
                'set(me,''checked'',''on'');',...
                'set(gcf,''UserData'',ud);',...
                'ppset = findobj(''name'',''matpplane Setup'');',...
                'sud = get(ppset,''UserData'');',...
                'sud.settings = ud.settings;',...
                'set(ppset,''UserData'',sud);',...
                'matpplane(''settings'');'];
            
            solver = dud.settings.solver;
            dpset.refine = 8;
            dpset.tol = dud.settings.tol;
            dpset.solver = 'Dormand Prince';
            dpset.stepsize = dud.settings.stepsize;
            dpset.hmax = 0;
            if strcmp(solver,'Dormand Prince')
                dpch = 'on';
            else
                dpch = 'off';
            end
            rk4set.refine = 1;
            rk4set.tol = dud.settings.tol;
            rk4set.solver = 'Runge-Kutta 4';
            rk4set.stepsize = dud.settings.stepsize;
            rk4set.hmax = 0;
            if strcmp(solver,'Runge-Kutta 4')
                rk4ch = 'on';
            else
                rk4ch = 'off';
            end
            ode15sset.refine = 1;
            ode15sset.tol = dud.settings.tol;
            ode15sset.solver = 'ode15s';
            ode15sset.stepsize = dud.settings.stepsize;
            ode15sset.hmax = 0;
            if strcmp(solver,'ode15s')
                ode15sch = 'on';
            else
                ode15sch = 'off';
            end
            ode23sset.refine = 1;
            ode23sset.tol = dud.settings.tol;
            ode23sset.solver = 'ode23s';
            ode23sset.stepsize = dud.settings.stepsize;
            ode23sset.hmax = 0;
            if strcmp(solver,'ode23s')
                ode23sch = 'on';
            else
                ode23sch = 'off';
            end
            ode113set.refine = 1;
            ode113set.tol = dud.settings.tol;
            ode113set.solver = 'ode113';
            ode113set.stepsize = dud.settings.stepsize;
            ode113set.hmax = 0;
            if strcmp(solver,'ode113')
                ode113ch = 'on';
            else
                ode113ch = 'off';
            end
            ode23set.refine = 1;
            ode23set.tol = dud.settings.tol;
            ode23set.solver = 'ode23';
            ode23set.stepsize = dud.settings.stepsize;
            ode23set.hmax = 0;
            if strcmp(solver,'ode23')
                ode23ch = 'on';
            else
                ode23ch = 'off';
            end
            ode45set.refine = 8;
            ode45set.tol = dud.settings.tol;
            ode45set.solver = 'ode45';
            ode45set.stepsize = dud.settings.stepsize;
            ode45set.hmax = 0;
            if strcmp(solver,'ode45')
                ode45ch = 'on';
            else
                ode45ch = 'off';
            end
            ode23tset.refine = 1;
            ode23tset.tol = dud.settings.tol;
            ode23tset.solver = 'ode23t';
            ode23tset.stepsize = dud.settings.stepsize;
            ode23tset.hmax = 0;
            if strcmp(solver,'ode23t')
                ode23tch = 'on';
            else
                ode23tch = 'off';
            end
            ode23tbset.refine = 1;
            ode23tbset.tol = dud.settings.tol;
            ode23tbset.solver = 'ode23tb';
            ode23tbset.stepsize = dud.settings.stepsize;
            ode23tbset.hmax = 0;
            if strcmp(solver,'ode23tb')
                ode23tbch = 'on';
            else
                ode23tbch = 'off';
            end
            
            dud.solver(1) = uimenu(mesolve,'label','Dormand Prince',...
                'checked',dpch,...
                'CallBack',solverstr,'UserData',dpset);
            
            dud.solver(2) = uimenu(mesolve,'label','Runge-Kutta 4',...
                'checked',rk4ch,...
                'CallBack',solverstr,'UserData',rk4set);
            
            dud.solver(3) = uimenu(mesolve,'label','ode45',...
                'checked',ode45ch,...
                'separator','on',...
                'CallBack',solverstr,'UserData',ode45set);
            
            dud.solver(4) = uimenu(mesolve,'label','ode23',...
                'checked',ode23ch,...
                'CallBack',solverstr,'UserData',ode23set);
            
            dud.solver(5) = uimenu(mesolve,'label','ode113',...
                'checked',ode113ch,...
                'CallBack',solverstr,'UserData',ode113set);
            
            dud.solver(6) = uimenu(mesolve,'label','ode15s',...
                'checked',ode15sch,...
                'CallBack',solverstr,'UserData',ode15sset);
            
            dud.solver(7) = uimenu(mesolve,'label','ode23s',...
                'checked',ode23sch,...
                'CallBack',solverstr,'UserData',ode23sset);
            
            dud.solver(8) = uimenu(mesolve,'label','ode23t',...
                'checked',ode23tch,...
                'CallBack',solverstr,'UserData',ode23tset);
            
            dud.solver(9) = uimenu(mesolve,'label','ode23tb',...
                'checked',ode23tbch,...
                'CallBack',solverstr,'UserData',ode23tbset);
            
            plotch = [
                'dud = get(gcf,''UserData'');',...
                'aud = get(dud.axes,''UserData'');',...
                'if aud.plot,',...
                '   aud.plot = 0;',...
                '   set(gcbo,''label'',''Plot while computing'');',...
                'else,',...
                '   aud.plot = 1;',...
                '   set(gcbo,''label'',''Do not plot while computing'');',...
                'end,',...
                'set(dud.axes,''UserData'',aud);'];
            
            medir = uimenu(menu(1),'label','Solution direction.');
            
            directionstr = ['ud = get(gcf,''UserData'');',...
                'me = gcbo;',...
                'ud.dir = get(me,''UserData'');',...
                'set(ud.direction,''checked'',''off'');',...
                'set(me,''checked'',''on'');',...
                'set(gcf,''UserData'',ud);'];
            
            dud.direction(1) = uimenu(medir,'label','Both',...
                'checked','on',...
                'UserData',0,...
                'CallBack',directionstr);
            dud.dir = 0;
            
            dud.direction(2) = uimenu(medir,'label','Forward',...
                'UserData',1,...
                'CallBack',directionstr);
            
            dud.direction(3) = uimenu(medir,'label','Back',...
                'UserData',-1,...
                'CallBack',directionstr);
            
            
            markstr = ['ud = get(gcf,''UserData'');',...
                'me = gcbo;',...
                'chkd = get(me,''checked'');',...
                'if strcmp(chkd,''on''),',...
                '  set(me,''checked'',''off'');',...
                '  ud.markflag = 0;',...
                'else,',...
                '  set(me,''checked'',''on'');',...
                '  ud.markflag = 1;',...
                'end,',...
                'set(gcf,''UserData'',ud);'];
            
            dud.markflag = 0;
            dud.mark = uimenu(menu(1),'label','Mark initial points.',...
                'checked','off',...
                'CallBack',markstr);
            
            [~] = uimenu(menu(1),'label','Export solution data.',...
                'CallBack','matpplane(''export'')',...
                'separator','off','tag','dexp');
            
            meplot = uimenu(menu(1),'label','Do not plot while computing',...
                'CallBack',plotch,'separator','on');
            
            menu(6) = uimenu(menu(1),'label','Make the Display Window inactive.',...
                'CallBack','matpplane(''hotcold'')','separator','on');
            
            dud.menu = menu;
            
            % View menu
            
            set(findobj(gcf,'label','&View'),'pos',6);
            set(findobj(gcf,'label','&Figure Toolbar'),...
                'CallBack','matpplane(''showbar'')');
            
            set(0,'showhiddenhandles',hhsetup);
            
            set(gcf,'WindowButtonDownFcn','matpplane(''down'')');
            set(ppdisp,'WindowButtonMotionFcn','matpplane(''cdisp'')');
            hh1 = [dud.axes,dud.title.axes,nframe,dud.notice,dbutt([1 2 3])];
            set(hh1,'units','norm');
            hh2 = [nframe,dud.notice,dbutt(2:3),dud.axes,dud.menu(1),meplot,megraph];
            set(hh2,'visible','on');
            
            set(ppdisp,'vis','on');
            dud.printstr = 'print -noui';
        end   % if (~isempty(ppdisp))  & else
        ppdispa = dud.axes;
        axes(ppdispa);
        cla
        xlabel(txname);
        ylabel(tyname);
        % Initialize the window matrix.
        
        set(findobj('tag','zbmenu'),'enable','off');
        tstr1 = [txname,' '' = ', txderstr];
        tstr2 = [tyname,' '' = ', tyderstr];
        tstr = char(tstr1,tstr2);
        dud.tstr = tstr;
        
        
        
        k = find(pflag);
        if ~isempty(k)
            lk = length(k);
            switch lk
                case 1
                    pstr1 = [pstring(k);{' '}];
                case 2
                    pstr1 = pstring(k);
                case 3
                    pstr1 = pstring(k([2,3]));
                    pstr2 = [pstring(k(1));{' '}];
                case 4
                    pstr1 = pstring(k([2,4]));
                    pstr2 = pstring(k([1,3]));
                case 5
                    pstr1 = pstring(k([3,5]));
                    pstr2 = pstring(k([2,4]));
                    pstr3 = [pstring(k(1));{' '}];
                case 6
                    pstr1 = pstring(k([3,6]));
                    pstr2 = pstring(k([2,5]));
                    pstr3 = pstring(k([1,4]));
            end
        end
        set(dud.title.eq,'string',tstr);
        set(dud.title.p1,'string',pstr1);
        ext = get(dud.title.p1,'extent');
        pos = get(dud.title.p1,'pos');
        p1 = min(.9, .93 - ext(3));
        pos(1) = p1;
        set(dud.title.p1,'pos',pos);
        set(dud.title.p2,'string',pstr2);
        ext = get(dud.title.p2,'extent');
        pos = get(dud.title.p2,'pos');
        p2 = min(.8, p1 - ext(3)-0.02);
        pos(1) = p2;
        set(dud.title.p2,'pos',pos);
        set(dud.title.p3,'string',pstr3);
        ext = get(dud.title.p3,'extent');
        pos = get(dud.title.p3,'pos');
        pos(1) = min(.7, p2 - ext(3)-0.02);
        set(dud.title.p3,'pos',pos);
        
        % Initialize important information as UserData data.
        
        dud.function = dfcn;
        dud.solhand = [];	% Handles to solution curves.
        dud.ephand = [];	% Handles to equilibrium points.
        dud.arr = [];	% Handles for the direction and vector fields.
        dud.eqpts = [];	% Equilibrium point data.
        dud.ics = [];      % Marked initial points.
        dud.wmat = [];
        dud.color = sud.color;
        set(ppdisp,'UserData',dud);
        if strcmp(dud.syst.fieldtype,'nullclines')
            menunull = findobj('tag','null');
            set(menunull,'label','Hide nullclines.');
        end
        ud.y = zeros(2,1);
        ud.i = 0;
        ud.line = 0;
        wind = dud.syst.wind(:);
        dwind = [wind(1); wind(3); -wind(2); -wind(4)];
        ud.DY = [wind(2)-wind(1); wind(4)-wind(3)];
        ud.cwind = dwind - dud.settings.magn*[ud.DY;ud.DY];
        ud.R = zeros(2,2);
        ud.rr = zeros(2,2);
        ud.perpeps = 0;
        ud.paraeps = 0;
        ud.sinkeps = 0;
        ud.minNsteps = 0;
        ud.turn = zeros(2,10);
        ud.tk = 0;
        ud.stop = 0;
        ud.gstop = 1;
        ud.plot = 1;
        set(dud.axes,'UserData',ud);
        tc = findobj('tag','testcase');
        if ~isempty(tc)
            if strcmp(dud.syst.name,'default system')
                set(tc,'vis','on');
            else
                set(tc,'vis','off');
            end
        end
        ppkbd = findobj('name','matpplane Keyboard input','vis','on');
        if ~isempty(ppkbd),matpplane('kbd'),end
        pplevel = findobj('name','matpplane Level sets');
        if ~isempty(pplevel),delete(pplevel),end
        set(dud.title.axes,'handle','off');
        matpplane('dirfield',ppdisp);
    end % if sud.flag == 1   &  else
    
elseif strcmp(action,'linear')
    
    % Initialize linearization window.
    
    % Linear takes the information from the EQPT window and initializes
    % the Linear Display window if it already exists.  If the Linear
    % Display window does not exist, it builds one.
    
    ppeqpt = gcf;
    
    % Get the information from matpplane Equilibrium ...
    
    sud = get(ppeqpt,'UserData');
    %WINvect = [-1 1 -1 1];
    jac = sud.jac;
    type = sud.type;
    %vectors = sud.vectors;
    ppdisp = findobj('name','matpplane Display');
    pdud = get(ppdisp,'UserData');
    settings = pdud.settings;
    system = sud.system;
    
    % Find matpplane Linearization if it exists.
    % If matpplane Linearization exists, update it.  If it does not build it.
    
    pplin= findobj('name','matpplane Linearization');
    if (~isempty(pplin))
        figure(pplin);
        dud = get(pplin,'UserData');
        dud.syst = system;
        dud.settings = settings;
        dfcn = dud.function;
    else
        pplin = figure('name','matpplane Linearization',...
            'numb','off',...
            'interrupt','on',...
            'visible','off',...
            'tag','matpplane');
        matpplane('figdefault',pplin);
        dud = get(pplin,'UserData');
        dud.dir = 0;
        dud.syst = system;
        dud.settings = settings;
        dfcn = '';
        fs = dud.fontsize;
        ssize = dud.ssize;
        r = ssize/10;
        ppaxw = 300;	% Default axes width
        ppaxh = 300;	% Default axes height
        ppaxl = 45;	% Default axes left
        %buttw = 40;     % Default button width
        titleh = 45;    % Default title height
        ppaxb = 35;
        ppdh = ppaxb + ppaxh + titleh;
        uni = get(0,'units');
        set(0,'units','pixels');
        ss = get(0,'screensize');
        set(0,'units',uni);
        sw = ss(3);sh = ss(4);
        if r*ppdh > sh -80
            r = (sh-80)/ppdh;
            fs = fs*r;
            lw = 0.5*r;
            set(gcf,'defaultaxesfontsize',fs,'defaultaxeslinewidth',lw);
            set(gcf,'defaulttextfontsize',fs);
            set(gcf,'defaultlinelinewidth',lw);
            set(gcf,'defaultuicontrolfontsize',fs*0.9);
        end
        
        axpos = r*[ppaxl,ppaxb,ppaxw,ppaxh];
        
        % Set up the plot axes.
        
        dud.axes = axes('units','pix',...
            'position',axpos,...
            'next','add',...
            'box','on',...
            'interrupt','on',...
            'xgrid','on',...
            'ygrid','on',...
            'drawmode','fast',...
            'visible','off',...
            'tag','linear axes');
        
        % Set up the title.
        
        dud.title.axes = axes('box','off','xlim',[0 1],'ylim',[0 1],...
            'units','pix','vis','off','xtick',-1,'ytick',-1,...
            'xticklabel','','yticklabel','');
        
        dud.title.eq = text(0.07,0.5,' ','vert','middle');
        dud.title.p1 = text(0.75,0.5,' ','vert','middle');
        dud.title.p2 = text(0.65,0.5,' ','vert','middle');
        tstr = {'x_2';'y^2'};
        set(dud.title.eq,'string',tstr,'units','pix');
        ext = get(dud.title.eq ,'extent');
        titleh = ext(4)+15*r;
        set(dud.title.eq,'units','data');
        taxpos = r*[ppaxl,ppaxb+ppaxh,ppaxw,titleh];
        set(dud.title.axes,'pos',taxpos,'color',get(gcf,'color'));
        
        % Finish the positions.
        
        % 		buttw = r*buttw;
        % 		butth = fs+10*r;
        % 		buttl = ppaxl + ppaxw + 5;
        ppdw = r*(ppaxl + ppaxw + 35);
        ppdh = r*(ppaxb+ppaxh)+titleh;
        
        ppdleft = sw-ppdw-60;
        ppdbot = 60;
        ppdpos = [ppdleft,ppdbot,ppdw,ppdh];
        set(pplin,'resize','on');
        set(pplin,'pos',ppdpos);
        
        %Arrflag = system.fieldtype;
        
        axpos = r*[ppaxl,ppaxb,ppaxw,ppaxh];
        
        stopstr = 'aud = get(gca,''UserData'');aud.stop = 4;set(gca,''UserData'',aud);';
        
        stoppos = [axpos(1)+axpos(3)-r*10, axpos(2)+axpos(4)+r*10, r*40,fs+10*r];
        
        dud.butt  = uicontrol('style','push',...
            'pos',stoppos,...
            'string','Stop','CallBack',stopstr,'vis','off','tag','stop');
        
        hhsetup = get(0,'showhiddenhandles');
        set(0,'showhiddenhandles','on');
        
        % Configure the Toolbar.
        
        set(pplin,'ToolBar','none');
        
        tmenu = findobj(gcf,'label','&Tools');
        delete(tmenu);
        
        % File menu
        
        fmenu = findobj(gcf,'label','&File');
        delete(findobj(fmenu,'label','&New Figure'));
        delete(findobj(fmenu,'label','&Open...'));
        delete(findobj(fmenu,'label','Pre&ferences...'));
        set(findobj(fmenu,'label','&Close'),'pos',1);
        set(findobj(fmenu,'label','&Save'),'pos',2,...
            'separator','off');
        set(findobj(fmenu,'label','Save &As...'),'pos',3);
        set(findobj(fmenu,'label','&Export...'),...
            'pos',4);
        delete(findobj(fmenu,'label','Pre&ferences...'));
        set(findobj(fmenu,'label','Pa&ge Setup...'),'pos',5);
        set(findobj(fmenu,'label','Print Set&up...'),'pos',6);
        set(findobj(fmenu,'label','Print Pre&view...'),'pos',7);
        set(findobj(fmenu,'label','&Print...'),'pos',8);
        
        [~] = uimenu(fmenu,'label','Quit matpplane',...
            'CallBack','matpplane(''quit'')','separator','on','pos',9);
        
        % View menu
        
        set(findobj(gcf,'label','&Figure Toolbar'),...
            'CallBack','matpplane(''showbar'')');
        
        % Solutions menu
        
        solmenu = uimenu('label','Solutions','pos',3);
        
        [~] = uimenu(solmenu,'label','Keyboard input.','CallBack',...
            'matpplane(''kbd'')','vis','on');
        
        [~]   = uimenu(solmenu,'label','Plot several solutions.',...
            'CallBack','matpplane(''several'')');
        
        
        fundCallBack = ['dud = get(gcf,''UserData'');',...
            'col = dud.color.sep;',...
            'vect = dud.vectors;',...
            'h = zeros(1,2);',...
            'h(1) = plot(2*[vect(1,1),-vect(1,1)],2*[vect(2,1),-vect(2,1)]);',...
            'h(2) = plot(2*[vect(1,2),-vect(1,2)],2*[vect(2,2),-vect(2,2)]);',...
            'set(h,''color'',col);',...
            'dud.solhand = [dud.solhand;h(:)];'...
            'set(gcf,''UserData'',dud);'];
        
        dud.menu = uimenu(solmenu,'CallBack',fundCallBack);
        
        markstr = ['ud = get(gcf,''UserData'');',...
            'me = gcbo;',...
            'chkd = get(me,''checked'');',...
            'if strcmp(chkd,''on''),',...
            '  set(me,''checked'',''off'');',...
            '  ud.markflag = 0;',...
            'else,',...
            '  set(me,''checked'',''on'');',...
            '  ud.markflag = 1;',...
            'end,',...
            'set(gcf,''UserData'',ud);'];
        
        dud.markflag = pdud.markflag;
        chkd = 'off';
        if dud.markflag
            chkd = 'on';
        end
        dud.mark = uimenu(solmenu,'label','Mark initial points.',...
            'checked',chkd,...
            'CallBack',markstr);
        % Edit menu
        
        emenu = findobj(gcf,'label','&Edit');
        
        [~] = uimenu(emenu,'label','Erase all solutions.',...
            'CallBack','matpplane(''dallsol'')',...
            'pos',1);
        
        [~] = uimenu(emenu,'label','Erase all marked initial points.',...
            'CallBack','matpplane(''dallics'')',...
            'separator','off',...
            'pos',2);
        
        [~] = uimenu(emenu,'label','Erase all graphics objects.',...
            'CallBack','matpplane(''dall'')',...
            'separator','off',...
            'pos',3);
        [~] = uimenu(emenu,'label','Delete a graphics object.',...
            'CallBack','matpplane(''delete'')',...
            'visible','on',...
            'pos',4);
        
        [~] = uimenu(emenu,'label','Enter text on the Display Window.',...
            'CallBack','matpplane(''text'')',...
            'pos',5);
        set(findobj(emenu,'label','&Undo'),'separator','on',...
            'pos',6);
        set(findobj(emenu,'label','Cu&t'),'pos',7);
        set(findobj(emenu,'label','&Copy'),'pos',8);
        set(findobj(emenu,'label','&Paste'),'pos',9);
        set(findobj(emenu,'label','Clea&r'),'pos',10);
        set(findobj(emenu,'label','&Select All'),'pos',11);
        set(findobj(emenu,'label','Copy &Figure'),'pos',12);
        set(findobj(emenu,'label','Copy &Options'),'pos',13);
        set(findobj(emenu,'label','F&igure Properties'),'pos',14);
        set(findobj(emenu,'label','&Axes Properties'),'pos',15);
        set(findobj(emenu,'label','C&urrent Object Properties'),'pos',16);
        
        % Options menu
        
        optmenu = uimenu('label','Options','visible','off');
        
        [~] = uimenu(optmenu,'label','Make the Display Window inactive.',...
            'CallBack','matpplane(''hotcold'')','separator','on');
        
        set(0,'showhiddenhandles',hhsetup);
        
        %       set(pplin,'WindowButtonDownFcn',['pplane9(''down'',',num2str(get(pplin,'Number')),')']);
        set(pplin,'WindowButtonDownFcn',['matpplane(''down'',',num2str(get(pplin,'Number')),')']);
        hh1 = [dud.axes,dud.title.axes];
        set(hh1,'units','norm');
        set(dud.axes,'visible','on');
        set(pplin,'vis','on');
        dud.printstr = 'print';
    end   % if (~isempty(pplin))  & else
    dud.ics = [];
    
    switch type
        case 1
            vstr = 'on';
            lstr = 'Plot stable and unstable orbits.';
        case {2,3}
            vstr = 'on';
            lstr = 'Plot fundamental modes.';
        otherwise
            vstr = 'off';
            lstr='';
    end
    set(dud.menu,'label',lstr,'vis',vstr);
    
    pplina = dud.axes;
    dud.vectors = sud.vectors;
    axes(pplina);
    cla
    xlabel('u');
    ylabel('v');
    
    % The title strings.
    
    tstr = {'u'' = A u + B v';'v'' = C u + D v'};
    dud.tstr = tstr;
    pstr2 = {['A = ',num2str(jac(1,1))];['C = ',num2str(jac(2,1))]};
    pstr1 = {['B = ',num2str(jac(1,2))];['D = ',num2str(jac(2,2))]};
    
    set(dud.title.eq,'string',tstr);
    set(dud.title.p1,'string',pstr1);
    ext = get(dud.title.p1,'extent');
    pos = get(dud.title.p1,'pos');
    p1 = min(.9, .93 - ext(3));
    pos(1) = p1;
    set(dud.title.p1,'pos',pos);
    set(dud.title.p2,'string',pstr2);
    ext = get(dud.title.p2,'extent');
    pos = get(dud.title.p2,'pos');
    pos(1) = min(.8, p1 - ext(3)-0.02);
    set(dud.title.p2,'pos',pos);
    
    if ~strcmp(dfcn,'') && exist(dfcn,'file')
        delete([tempdir,dfcn,'.m']);
    end
    tee = clock;
    tee = ceil(tee(6)*100);
    dfcn=['pptp',num2str(tee)];
    fcnstr = ['function ypr = ',dfcn,'(t,y)\n\n'];
    commstr = '%%%% Created by matpplane\n\n';
    astr = ['A = [',num2str(jac(1,1)),',',num2str(jac(1,2)),';',...
        num2str(jac(2,1)),',',num2str(jac(2,2)),'];\n'];
    dstr = 'ypr = A * y;';
    ppf = fopen([tempdir,dfcn,'.m'],'w');
    fprintf(ppf,fcnstr);
    fprintf(ppf,commstr);
    fprintf(ppf,astr);
    fprintf(ppf,dstr);
    fclose(ppf);
    
    % Initialize important information as UserData data.
    
    dud.function = dfcn;
    dud.solhand = [];	% Handles to solution curves.
    dud.ephand = [];	% Handles to equilibrium points.
    dud.contours = [];
    dud.arr = [];		% Handles for the direction and vector
    % fields.
    dud.eqpts = [];		% Equilibrium point data.
    dud.notice = 0;
    dud.syst.wind = [-1 1 -1 1]';
    dud.color = sud.color;
    ud.y = zeros(2,1);
    ud.i = 0;
    ud.line = 0;
    wind = [-1 1 -1 1]';
    dwind = [wind(1); wind(3); -wind(2); -wind(4)];
    ud.DY = [wind(2)-wind(1); wind(4)-wind(3)];
    ud.cwind = dwind - dud.settings.magn*[ud.DY;ud.DY];
    ud.R = zeros(2,2);
    ud.rr = zeros(2,2);
    ud.perpeps = 0;
    ud.paraeps = 0;
    ud.sinkeps = 0;
    ud.turn = zeros(2,10);
    ud.tk = 0;
    ud.stop = 0;
    ud.gstop = 1;
    ud.plot = 1;
    dud.ephand = plot(0,0,...
        'color',dud.color.eqpt,...
        'markersize',20,...
        'marker','.');
    set(dud.axes,'UserData',ud);
    set(pplin,'UserData',dud);
    ppkbd = findobj('name','matpplane Keyboard input','vis','on');
    if ~isempty(ppkbd),matpplane('kbd'),end
    matpplane('dirfield',pplin);
    
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 matpplane Display and get the UserData data.
    
    disph = input1;   % This could be ppdisp or pplin.
    
    dud = get(disph,'UserData');
    color = dud.color;
    
    dfcn = dud.function;
    ppdispa = dud.axes;
    WINvect = dud.syst.wind;
    %settings = dud.settings;
    notice = dud.notice;
    
    
    % NPK attempt to setup the notice window for linearization:
    %    nframe = uicontrol('style','frame','visible','on');
    nstr = {'More';'than';'five';'lines';'of text'};
    dud.notice = uicontrol('style','text',...
        'horiz','left',...
        'string',nstr,'visible','on');
    dud.noticeflag = 1;
    
    
    if ishghandle(notice)
        
        try
            nstr = get(notice,'string');
            nstr(1:4)=nstr(2:5);
            nstr{5,1} = 'Computing the field elements.';
            set(notice,'string',nstr);
        catch
        end
        
        
        % Augment the window matrix
        dud.wmat = [];
        wmat = dud.wmat;
        wrows = size(wmat,1);
        wflag = 0;
        for k = 1:wrows
            if wmat(k,:)==WINvect
                wflag = 1;
            end
        end
        if wflag == 0
            wmat = [wmat;WINvect];
            dud.wmat = wmat;
        end
        if wrows
            hhsetup = get(0,'showhiddenhandles');
            set(0,'showhiddenhandles','on');
            hhh = findobj('tag','zbmenu');
            set(hhh,'enable','on');
            set(0,'showhiddenhandles',hhsetup);
        end
    end
    
    
    Xmin = WINvect(1);
    Xmax = WINvect(2);
    Ymin = WINvect(3);
    Ymax = WINvect(4);
    
    N = dud.syst.npts;
    if strcmp(get(disph,'name'),'matpplane Linearization')
        N = floor(3*N/4);
    end
    deltax=(Xmax - Xmin)/(N-1);
    deltay=(Ymax - Ymin)/(N-1);
    
    % Set up the display window.
    
    Dxint=[Xmin-deltax,Xmax+deltax];
    Dyint=[Ymin-deltay,Ymax+deltay];
    
    % Set up the original mesh.
    
    XXXg=Xmin + deltax*(0:N-1);
    YYYg=Ymin + deltay*(0:N-1);
    
    [Xx,Yy]=meshgrid(XXXg,YYYg);
    
    % Calculate the line and vector fields.
    
    Xx=Xx(:);Yy=Yy(:);
    %Ww = zeros(size(Xx)); % Unused preallocated variable
    Ww = feval(dfcn,0,[Xx';Yy']);
    Vv = Ww(1,:) + Ww(2,:)*1i;
    Vv = Vv.';
    Arrflag = dud.syst.fieldtype;
    
    
    mgrid = Xx+Yy.*1i; % 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.');
    [~,c] = size(aa1);
    aa1=[aa1;NaN*ones(1,c)];
    aa1=aa1(:);
    
    arrow = [0,1,.7,1,.7].' + [0,0,.25,0,-.25].' * 1i;
    zz=sign(zz).*((abs(zz)).^(1/3));
    scale = 0.9*sc./max(max(abs(zz)));
    aa2 = scale*arrow*zz +ones(size(arrow))*(mgrid.');
    [~,c] = size(aa2);
    aa2=[aa2;NaN*ones(1,c)];
    aa2=aa2(:);
    axes(ppdispa);
    
    arr = dud.arr;	% Delete the old field data.
    if isstruct(arr)
        hand = [arr.hx;arr.hy;arr.lines;arr.arrows;arr.barrows];
        delete(hand);
    end
    NN = N;
    N = 50; k = ceil(N/NN);
    deltax=(Xmax - Xmin)/(N-1);
    deltay=(Ymax - Ymin)/(N-1);
    
    % Set up the original mesh.
    
    XXXg=Xmin + deltax*(-k:N+k);
    YYYg=Ymin + deltay*(-k:N+k);
    
    [Xx,Yy]=meshgrid(XXXg,YYYg);
    
    % Calculate the line and vector fields.
    
    Xxx=Xx(:);Yyy=Yy(:);
    %Ww = zeros(size(Xxx)); %Unused preallocated variable
    Ww = feval(dfcn,0,[Xxx';Yyy']);
    DX = Ww(1,:)';
    DY = Ww(2,:)';
    minx = min(DX);
    maxx = max(DX);
    miny = min(DY);
    maxy = max(DY);
    DR = Xxx.*DX + Yyy.*DY;
    DTheta = Xxx.*DY - Yyy.*DX;
    minr = min(DR);
    maxr = max(DR);
    minth = min(DTheta);
    maxth = max(DTheta);
    DX = reshape(DX,N+2*k+1,N+2*k+1);
    DY = reshape(DY,N+2*k+1,N+2*k+1);
    DR = reshape(DR,N+2*k+1,N+2*k+1);
    DTheta = reshape(DTheta,N+2*k+1,N+2*k+1);
    if minx < 0 && 0 < maxx
        [~,hx] = contour(Xx,Yy,DX,[0,0],'--');
        set(hx,'visible','off','color',color.xcline,'linestyle','--');
    else
        hx = zeros(0,1);
    end
    if miny < 0 && 0 < maxy
        [~,hy] = contour(Xx,Yy,DY,[0,0],'--');
        set(hy,'visible','off','color',color.ycline,'linestyle','--');
    else
        hy = zeros(0,1);
    end
    if minr < 0 && 0 < maxr
        [~,hr] = contour(Xx,Yy,DR,[0,0],'--');
        set(hr,'visible','off','color',color.xcline,'linestyle','--');
    else
        hr = zeros(0,1);
    end
    if minth < 0 && 0 < maxth
        [~,hth] = contour(Xx,Yy,DTheta,[0,0],'--');
        set(hth,'visible','off','color',color.ycline,'linestyle','--');
    else
        hth = zeros(0,1);
    end
    arrh1 = plot(real(aa1),imag(aa1),'color',color.arrows,'visible','off');
    arrh2 = plot(real(aa2),imag(aa2),'color',color.arrows,'visible','off');
    
    % We plot both the line field and the vector field.  Then we
    % control which is seen by manipulating the visibility.
    Xmin = Xmin - k*deltax;
    Xmax = Xmax + k*deltax;
    Ymin = Ymin - k*deltay;
    Ymax = Ymax + k*deltay;
    %     Zz = (Xx-Xmin).*(Xmax-Xx).*(Yy-Ymin).*(Ymax-Yy).*abs(DX).*abs(DY);
    %     Zzc = Zz(2:N+k+2,2:N+k+2);
    %     Xxc = Xx(2:N+k+2,2:N+k+2);
    %     Yyc = Yy(2:N+k+2,2:N+k+2);
    %     DXc = DX(2:N+k+2,2:N+k+2);
    %     DYc = DY(2:N+k+2,2:N+k+2);
    %     Zzl = Zz(1:N+k+1,2:N+k+2);
    %     Zzr = Zz(3:N+k+3,2:N+k+2);
    %     Zzd = Zz(2:N+k+2,1:N+k+1);
    %     Zzu = Zz(2:N+k+2,3:N+k+3);
    %     Zzdl = Zz(1:N+k+1,1:N+k+1);
    %     Zzdr = Zz(3:N+k+3,1:N+k+1);
    %     Zzul = Zz(1:N+k+1,3:N+k+3);
    %     Zzur = Zz(3:N+k+3,3:N+k+3);
    %kk = find((Zzc>Zzl)&(Zzc>Zzr)&(Zzc>Zzd)&(Zzc>Zzu)...
    %    &(Zzc>Zzdl)&(Zzc>Zzdr)&(Zzc>Zzul)&(Zzc>Zzur));
    %Xxx = Xxc(kk);
    %Yyy = Yyc(kk);
    %DXx = DXc(kk);
    %DYy = DYc(kk);
    ll = length(hx);
    %Xxx = Xxx(:);
    %Yyy = Yyy(:);
    %DXx = DXx(:);
    %DYy = DYy(:);
    %if ~exist('both')   %#ok<EXIST> % I don't think 'both' ever exists so this always gets called
    Xxx = zeros(0,1);
    Yyy = zeros(0,1);
    DXx = zeros(0,1);
    DYy = zeros(0,1);
    %end
    for j = 1:ll
        xd = get(hx(j),'xdata');xd = xd(:);
        yd = get(hx(j),'ydata');yd = yd(:);
        zxd = feval(dfcn,0,[xd';yd']);
        yxd = zxd(2,:);yxd = yxd(:);
        ayxd = abs(yxd).*(xd-Xmin).*(Xmax-xd).*(yd-Ymin).*(Ymax-yd);
        NNN = length(ayxd);
        ayxdc = ayxd(2:NNN-1);
        yxdc = yxd(2:NNN-1);
        ayxdl = ayxd(1:NNN-2);
        ayxdr = ayxd(3:NNN);
        lll = find((ayxdc>ayxdl)&(ayxdc>ayxdr));
        Xx = xd(lll+1);
        Yy = yd(lll+1);
        Xxx = [Xxx;Xx(:)]; %#ok<AGROW>
        Yyy = [Yyy;Yy(:)]; %#ok<AGROW>
        yxdcp = yxdc(lll);
        yxdcp = yxdcp(:);
        DXx = [DXx;zeros(size(yxdcp))]; %#ok<AGROW>
        DYy = [DYy;yxdcp]; %#ok<AGROW>
    end
    ll = length(hy);
    for j = 1:ll
        xd = get(hy(j),'xdata');xd = xd(:);
        yd = get(hy(j),'ydata');yd = yd(:);
        zxd = feval(dfcn,0,[xd';yd']);
        yxd = zxd(1,:);yxd = yxd(:);
        ayxd = abs(yxd).*(xd-Xmin).*(Xmax-xd).*(yd-Ymin).*(Ymax-yd);
        NNN = length(ayxd);
        ayxdc = ayxd(2:NNN-1);
        yxdc = yxd(2:NNN-1);
        ayxdl = ayxd(1:NNN-2);
        ayxdr = ayxd(3:NNN);
        lll = find((ayxdc>ayxdl)&(ayxdc>ayxdr));
        Xx = xd(lll+1);
        Yy = yd(lll+1);
        Xxx = [Xxx;Xx(:)]; %#ok<AGROW>
        Yyy = [Yyy;Yy(:)]; %#ok<AGROW>
        yxdcp = yxdc(lll);
        yxdcp = yxdcp(:);
        DXx = [DXx;yxdcp]; %#ok<AGROW>
        DYy = [DYy;zeros(size(yxdcp))]; %#ok<AGROW>
    end
    mgrid = Xxx(:) + Yyy(:)*1i;
    Vv = DXx(:) + DYy(:)*1i;
    Vv = sign(Vv);       %.*((abs(Vv)).^(1/3));
    sc = (N/20)*min(deltax,deltay);
    aa3 = sc*arrow*(Vv.') +ones(size(arrow))*(mgrid.');
    [~,c] = size(aa3);
    aa3 = [aa3;NaN*ones(1,c)];
    aa3=aa3(:);
    plot(real(aa3),imag(aa3),'color',color.narrows,'visible','off');
    
    switch Arrflag
        case 'nullclines'
            %      set([hx;hy;arrh3],'vis','on');
            set([hx;hy],'vis','on');
        case 'lines'
            set(arrh1,'visible','on');
        case 'arrows'
            set(arrh2,'visible','on');
    end
    dud.arr.lines = arrh1;	   % Save the handles for later use.
    dud.arr.arrows = arrh2;	   % Save the handles for later use.
    dud.arr.hx = hx;
    dud.arr.hy = hy;
    %   dud.arr.barrows = arrh3;
    dud.arr.barrows = [];
    dud.arr.hr = hr;
    dud.arr.hth = hth;
    
    menull = findobj('tag','null');
    mernull = findobj('tag','rnull');
    if strcmp(get(menull,'label'),'Delete nullclines.')
        set([hx;hy],'vis','on');
    elseif strcmp(get(mernull,'label'),'Delete polar nullclines.')
        set([hr;hth],'vis','on');
    end
    try
        if ishghandle(notice)
            nstr = get(notice,'string');
            nstr(1:4) = nstr(2:5);
            nstr{5,1} = 'Ready.';
            set(notice,'string',nstr);
        end
    catch
    end
    set(disph,'UserData',dud);
    axis([Dxint,Dyint]);
    
elseif strcmp(action,'hotcold')
    
    % 'hotcold' is the CallBackback for the menu selection that makes the
    % Display Window active or inactive.
    
    ppdisp = gcf;
    dud = get(ppdisp,'UserData');
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    mehc = dud.menu(6);
    mehclabel=get(mehc,'label');
    if contains(mehclabel,'inactive')
        set(ppdisp,'WindowButtonDownFcn',' ');
        set(mehc,'label','Make the Display Window active.');
        nstr{5,1} = 'The Display Window is not active.';
        set(dud.notice,'string',nstr);
    else
        set(ppdisp,'WindowButtonDownFcn','matpplane(''down'')');
        set(mehc,'label','Make the Display Window inactive.');
        nstr{5,1} = 'The Display Window is active.';
        set(dud.notice,'string',nstr);
    end
    
elseif strcmp(action,'down')
    
    % 'down' is the Window Button Down CallBack.  It starts the computation of
    % solutions from a click of the mouse.
    
    disph = gcf;
    seltype = get(disph,'selectiontype');
    if strcmp(seltype,'alt')
        matpplane('zoom');
        return
    elseif strcmp(seltype,'extend')
        matpplane('zoomsqd');
        return
    end
    dud = get(disph,'UserData');
    ax = dud.axes;
    ch = findobj('type','uicontrol','enable','on');
    set(ch,'enable','inactive');
    wbdf = get(disph,'WindowbuttonDownFcn');
    set(disph,'WindowbuttonDownFcn','');
    axes(ax);
    initpt = get(ax,'currentpoint');
    initpt = initpt(1,[1,2]);
    if dud.markflag
        h = plot(initpt(1),initpt(2),'.k');
        dud.ics = [dud.ics,h];
        set(disph,'UserData',dud);
    end
    matpplane('solution',initpt,disph);
    set(disph,'WindowbuttonDownFcn',wbdf);
    %  set([ch;mh],'enable','on');
    set(ch,'enable','on');
    notice = dud.notice;
    if ishghandle(notice)
        nstr = get(notice,'string');
        nstr(1:4) = nstr(2:5);
        nstr{5,1} = 'Ready.';
        set(notice,'string',nstr)
    end
    
elseif strcmp(action,'several')
    
    % 'several' allows the UserData to pick several initial points at once.
    % This is not needed in X-windows, but it is on the Macintosh.
    
    disph = gcf;
    ch = findobj('type','uicontrol','enable','on');
    %  mh = findobj('type','uimenu','enable','on');
    set(ch,'enable','inactive');
    %  set(mh,'enable','off')
    wbdf = get(disph,'WindowbuttonDownFcn');
    set(disph,'WindowbuttonDownFcn','');
    dud = get(disph,'UserData');
    notice = dud.notice;
    if ishghandle(notice)
        nstr = get(notice,'string');
        nstr(1:4) = nstr(2:5);
        nstr{5,1} = 'Pick initial points with the mouse. Enter "Return" when finished.';
        set(notice,'string',nstr)
    end
    %   [X,Y]=ppginput;
    [X,Y]=ginput;
    NN = length(X);
    for k = 1:NN
        initpt = [X(k),Y(k)];
        matpplane('solution',initpt,disph);
    end
    if ishghandle(notice)
        nstr = get(notice,'string');
        nstr(1:4) = nstr(2:5);
        nstr{5,1} = 'Ready.';
        set(notice,'string',nstr);
    end
    set(ch,'enable','on');
    %  set([ch;mh],'enable','on');
    set(disph,'WindowbuttonDownFcn',wbdf);
    
elseif strcmp(action,'test case')
    
    tic
    ppdisp = gcf;
    for k = 1:10
        initpt = k*[-.2 .2];
        matpplane('solution',initpt,ppdisp);
        initpt = k*[-.2 -.2];
        matpplane('solution',initpt,ppdisp);
    end
    matpplane('solution',[1,-1],ppdisp);
    dud = get(gcf,'UserData');
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5,1} = 'Ready.';
    set(dud.notice,'string',nstr)
    toc
    
elseif strcmp(action,'solution')
    
    % 'solution' effects the computation and (erasemode == none) plotting of
    % solutions.  It also stores the data as appropriate.
    
    disph = input2;
    dud = get(disph,'UserData');
    %tcol = dud.color.temp;
    pcol = dud.color.orb;
    notice = dud.notice;
    initpt = input1(:);
    dfcn = dud.function;
    ppdispa = dud.axes;
    settings = dud.settings;
    ptstr = [' (',num2str(initpt(1),2), ', ', num2str(initpt(2),2), ')'];
    refine = settings.refine;
    %ssize = settings.stepsize;
    tol = settings.tol;
    ud = get(dud.axes,'UserData');
    rtol = tol;
    atol = tol*ud.DY*1e-4';
    
    
    if length(initpt)  == 2
        AA = -1e6;
        BB = 1e6;
        switch dud.dir
            case 0
                intplus = [0, BB];
                intminus = [0, AA];
            case -1
                intplus = [0, 0];
                intminus = [0, AA];
            case 1
                intplus = [0, BB];
                intminus = [0, 0];
        end
        
    else
        intplus = [initpt(3),initpt(5)];
        intminus = [initpt(3),initpt(4)];
        initpt = initpt(1:2);
        
    end
    stopbutt = findobj(disph,'tag','stop');
    
    set(stopbutt,'vis','on','enable','on');
    
    solver = settings.solver;
    switch solver
        case 'Dormand Prince'
            solh = @matppdp45;
            opt = disph;
        case 'Runge-Kutta 4'
            solh = @matpprk4;
            opt = disph;
        case 'ode45'
            solh = @ode45;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',rtol,'Abstol',atol);
        case 'ode23'
            solh = @ode23;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',rtol,'Abstol',atol);
        case 'ode113'
            solh = @ode113;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',rtol,'Abstol',atol);
        case 'ode15s'
            solh = @ode15s;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',rtol,'Abstol',atol);
        case 'ode23s'
            solh = @ode23s;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',rtol,'Abstol',atol);
        case 'ode23t'
            solh = @ode23t;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',rtol,'Abstol',atol);
        case 'ode23tb'
            solh = @ode23tb;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',rtol,'Abstol',atol);
    end
    
    %exist(dfcn);  % does nothing
    dfh = str2func(dfcn);
    cflag = 0;
    
    if intplus(2)>intplus(1)
        cflag = cflag + 1;
        if ishghandle(notice)
            nstr = get(notice,'string');
            nstr(1:4) = nstr(2:5);
            nstr{5} = ['The forward orbit from',ptstr];
            set(notice,'string',nstr);
        end
        drawnow
        
        %        [tp,xp] = feval(solh,dfh,intplus,initpt,opt);
        [tp,xp] = solh(dfh,intplus,initpt,opt);
        aud = get(ppdispa,'UserData');
        hnew1 = aud.line;
    end
    
    if intminus(2) < intminus(1)
        cflag = cflag + 2;
        if ishghandle(notice)
            nstr = get(notice,'string');
            nstr(1:4) = nstr(2:5);
            nstr{5} = ['The backward orbit from',ptstr];
            set(notice,'string',nstr);
        end
        drawnow
        %        [tm,xm] = feval(solh,dfh,intminus,initpt,opt);
        [tm,xm] = solh(dfh,intminus,initpt,opt);
        aud = get(ppdispa,'UserData');
        hnew2 = aud.line;
        
        set(stopbutt,'vis','off');
    end	  % if intminus(2) < intminus(1)
    
    
    
    % Store the trajectory.
    
    switch cflag
        case 1 % positive only
            set(hnew1,'xdata',xp(:,1),'ydata',xp(:,2),'zdata',tp,'color',pcol);
            set(hnew1,'erase','normal');
            dud.solhand = [dud.solhand;hnew1];
        case 2 % negative only
            x = flipud(xm);
            t = flipud(tm);
            set(hnew2,'xdata',x(:,1),'ydata',x(:,2),'zdata',t,'color',pcol);
            set(hnew2,'erase','normal');
            dud.solhand = [dud.solhand;hnew2];
            
        case 3 % both directions
            x = flipud(xm);
            t = flipud(tm);
            x=[x;xp];
            t=[t;tp];
            delete(hnew1);
            set(hnew2,'xdata',x(:,1),'ydata',x(:,2),'zdata',t,'color',pcol);
            set(hnew2,'erase','normal');
            dud.solhand = [dud.solhand;hnew2];
    end	 % switch cflag
    set(disph,'UserData',dud);
    
elseif strcmp(action,'kcompute')
    
    % 'kcompute' is the CallBack back for the Compute
    % button on the matpplane Keyboard figure.
    
    compute = 1;
    kh = get(gcf,'UserData');
    ppdisp = kh.fig;
    if (isempty(ppdisp))
        matpplane('confused');
    end
    dud = get(ppdisp,'UserData');
    ppdispa = dud.axes;
    aud = get(ppdispa,'UserData');
    ppset = findobj('name','matpplane Setup');
    sud = get(ppset,'UserData');
    ch = findobj('type','uicontrol','enable','on');
    set(ch,'enable','inactive');
    set(ppdisp,'WindowbuttonDownFcn','');
    xstr = get(kh.xval,'string');
    ystr = get(kh.yval,'string');
    pnameh = sud.h.pname;
    pvalh = sud.h.pval;
    %pflag = zeros(1,4);
    %perr = [];
    for kk = 1:6
        pn = char(get(pnameh(kk),'string'));
        pv = char(get(pvalh(kk),'string'));
        if ~isempty(pn)
            if isempty(pv)
                %perr = pvalh(kk);
            else
                xstr = matpplane('paraeval',pn,pv,xstr);
                ystr = matpplane('paraeval',pn,pv,ystr);
            end
        end
    end
    xvalue = str2double(xstr);
    yvalue = str2double(ystr);
    
    
    if get(kh.spec,'value')
        tzero = str2double(get(kh.tval,'string'));
        t0 = str2double(get(kh.t0,'string'));
        tf = str2double(get(kh.tf,'string'));
        initpt = [xvalue,yvalue,tzero,t0,tf];
        str1 = 'Values must be entered for all of the entries.';
        if (length(initpt) ~= 5)
            warndlg({str1},'Illegal input');
            compute = 0;
        elseif tf <= t0
            warndlg({'The final time of the computation interval';...
                'must be smaller than the initial time.'},'Illegal input');
            compute = 0;
        elseif (tzero < t0) || (tzero > tf)
            str2 = 'The initial time must be in the computation interval.';
            warndlg(str2,'Illegal input');
            compute = 0;
        end
        aud.gstop = 0;
        set(ppdispa,'UserData',aud);
    else
        initpt = [xvalue,yvalue];
        str1 = 'Values must be entered for all of the entries.';
        if (length(initpt) ~= 2)
            warndlg({str1},'Illegal input');
            compute = 0;
        end
    end	% if get(kh.spec,'value')
    
    if compute
        if dud.markflag
            figure(ppdisp)
            h = plot(initpt(1),initpt(2),'.k');
            dud.ics = [dud.ics,h];
            set(ppdisp,'UserData',dud);
        end
        matpplane('solution',initpt,ppdisp);
    end
    if dud.notice ~= 0
        nstr = get(dud.notice,'string');
        nstr(1:4) = nstr(2:5);
        nstr{5} = 'Ready.';
        set(dud.notice,'string',nstr);
    end
    %  set([ch;mh],'enable','on');
    set(ch,'enable','on');
    set(ppdisp,'WindowbuttonDownFcn','matpplane(''down'')');
    aud.gstop = 1;
    set(ppdispa,'UserData',aud);
    
    
elseif strcmp(action,'kbd')
    
    % 'kbd' is the CallBackback for the Keyboard Input menu selection.  It
    % sets up the matpplane Keyboard figure which allows accurate input of
    % initial values using the keyboard.
    
    ppdisp = gcf;  % The figure to be plotted in.
    dud = get(ppdisp,'UserData');
    Xname = dud.syst.xvar;
    Yname = dud.syst.yvar;
    xnstr = ['The initial value of ',Xname,' = '];
    ynstr = ['The initial value of ',Yname,' = '];
    tnstr = 'The initial value of t = ';
    ppkbd = findobj('name','matpplane Keyboard input');
    if  ~isempty(ppkbd)
        delete(ppkbd);
    end
    
    ppkbd = figure('name','matpplane Keyboard input',...
        'vis','off',...
        'numb','off','tag','matpplane');
    
    matpplane('figdefault',ppkbd);
    set(ppkbd,'menubar','none');
    
    kbd.fr1 = uicontrol('style','frame');
    
    kbd.fr2 = uicontrol('style','frame');
    
    kbd.fr3 = uicontrol('style','frame');
    
    dname = get(ppdisp,'name');
    kbd.which = uicontrol('style','text','horiz','center',...
        'string',['Data for ',dname]);
    
    kbd.inst = uicontrol('style','text','horiz','left',...
        'string','Enter the initial conditions:');
    
    kbd.xname = uicontrol('style','text',...
        'horiz','right','string',xnstr);
    
    kbd.xval = uicontrol('style','edit',...
        'string','','CallBack','');
    
    kbd.yname = uicontrol('style','text',...
        'horiz','right',...
        'string',ynstr);
    
    kbd.yval = uicontrol('style','edit',...
        'string','');
    
    kbd.tname = uicontrol('style','text',...
        'horiz','right',...
        'string',tnstr);
    
    kbd.tval = uicontrol('style','edit');
    
    kbd.t0 = uicontrol('style','edit');
    
    kbd.tf = uicontrol('style','edit');
    
    kbd.t = uicontrol('style','text','string','<= t<= ');
    
    kbd.spec = uicontrol('style','check','horiz','center',...
        'string','Specify a computation interval.');
    
    kbd.comp = uicontrol('style','push',...
        'string','Compute','CallBack','matpplane(''kcompute'')');
    
    kbd.close = uicontrol('style','push',...
        'string','Close','CallBack','set(gcf,''vis'',''off'')');
    kbd.fig = ppdisp;
    left = 5; varl = 70; % buttw = 60;
    frsep = 1;
    nudge = 3;
    xex = get(kbd.xname,'extent');
    ht = xex(4)+nudge;
    yex = get(kbd.yname,'extent');
    nwdth = max(xex(3),yex(3)) + nudge;
    varl = varl*ht/19;
    fr1bot = 2*left + ht;
    fr1ht = 4*nudge + 3*ht;
    frw = 2*nudge + nwdth + varl;
    fr1wind = [left,fr1bot,frw,fr1ht];
    set(kbd.fr1,'pos',fr1wind);
    tnb = fr1bot + nudge;
    nl = left+nudge;
    vl = nl + nwdth;
    tvwind = [vl,tnb,varl,ht];
    tnwind = [nl,tnb,nwdth,ht];
    set(kbd.tname,'pos',tnwind);
    set(kbd.tval,'pos',tvwind);
    intext = get(kbd.t,'extent');
    tw = intext(3);
    margin = (frw -tw - 2*varl)/2;
    t0l = left+margin;
    intbot = tnb +ht + nudge;
    t0wind = [t0l,intbot,varl,ht];
    set(kbd.t0,'pos',t0wind);
    tl = t0l + varl;
    tfl = tl + tw;
    twind = [tl,intbot,tw,ht];
    tfwind = [tfl,intbot,varl,ht];
    set(kbd.t,'pos',twind);
    set(kbd.tf,'pos',tfwind);
    specb = intbot + ht + nudge;
    specw = frw - 2*nudge;
    specwind = [nl,specb,specw,ht];
    set(kbd.spec,'pos',specwind);
    fr2bot = fr1bot + fr1ht +frsep;
    ynb = fr2bot + nudge;
    xnb = ynb + ht;
    xnwind = [nl,xnb,nwdth,ht];
    ynwind = [nl,ynb,nwdth,ht];
    xvwind = [vl,xnb,varl,ht];
    yvwind = [vl,ynb,varl,ht];
    instb = xnb + ht + nudge;
    instw = nwdth + varl;
    instwind = [nl,instb,instw,ht];
    fr2ht = 4*nudge + 3*ht;
    fr3bot = fr2bot + fr2ht+frsep;
    fr3ht = 2*nudge + ht;
    frw = 2*nudge + nwdth + varl;
    whichbot = fr3bot + nudge;
    whichw = frw - 2*nudge;
    whichwind = [nl,whichbot,whichw,ht];
    fr3wind = [left,fr3bot,frw,fr3ht];
    fr2wind = [left,fr2bot,frw,fr2ht];
    figw = 2*left + frw;
    fight = 3*left + ht + fr1ht + fr2ht + fr3ht;
    figwind = [30,300,figw,fight];
    buttw = (frw-left)/2;
    closel = left;
    compl = 2*left+buttw;
    clwind = [closel,left,buttw,ht];
    compwind = [compl,left,buttw,ht];
    set(ppkbd,'pos',figwind);
    set(kbd.fr2,'pos',fr2wind);
    set(kbd.fr3,'pos',fr3wind);
    set(kbd.which,'pos',whichwind);
    set(kbd.inst,'pos',instwind);
    set(kbd.xname,'pos',xnwind);
    set(kbd.yname,'pos',ynwind);
    set(kbd.xval,'pos',xvwind);
    set(kbd.yval,'pos',yvwind);
    set(kbd.comp,'pos',compwind);
    set(kbd.close,'pos',clwind);
    specCallBack = [
        'ud = get(gcf,''UserData'');',...
        'if get(gcbo,''value''),',...
        '  set([ud.t0,ud.t,ud.tf,ud.tname,ud.tval],''enable'',''on'');',...
        'else,',...
        '  set([ud.t0,ud.t,ud.tf,ud.tname,ud.tval],''enable'',''off'');',...
        'end'];
    
    set(kbd.spec,'CallBack',specCallBack);
    set(ppkbd,'resize','on');
    set(findobj(ppkbd,'type','uicontrol'),'units','normal');
    
    set([kbd.tval,kbd.t0],'string','0');
    set(kbd.spec,'value',0);
    set(ppkbd,'UserData',kbd,'vis','on');
    set([kbd.t0,kbd.t,kbd.tf,kbd.tname,kbd.tval],'enable','off')
    set(findobj(ppkbd,'type','uicontrol'),'units','normal');
    edith = findobj(ppkbd,'style','edit');
    set(edith,'backgroundcolor','w');
    figure(ppkbd)
    
elseif strcmp(action,'eqpt')
    
    % Find and classify equilibrium points.
    
    ppdisp = findobj('name','matpplane Display');
    dud = get(ppdisp,'UserData');
    col = dud.color.eqpt;
    dfcn = dud.function;
    %dbutt = dud.butt;
    %menu = dud.menu;
    ppdispa =dud.axes;
    Dx = get(ppdispa,'xlim');
    Dy = get(ppdispa,'ylim');
    DY = [Dx(2)-Dx(1);Dy(2)-Dy(1)];
    
    epsilon = 1e-4;
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'Choose an approximation with the mouse.';
    set(dud.notice,'string',nstr);
    
    %  z0 = ppginput(1);
    z0 = ginput(1);
    
    z = matpplane('newton',z0,dfcn);
    
    flag = z(:,4);
    A = z(:,2:3);
    B = real(A);
    k = find(abs(B)<1e-8);
    B(k) = zeros(size(k));
    C = imag(A);
    k = find(abs(C)<1e-8);
    C(k) = zeros(size(k));
    A = B+C*1i;
    z = z(:,1);
    
    if ~flag | norm((z-z0')./DY) > 1/5   %#ok<OR2> % flag not scalar, can't use ||
        nstr(1:4) = nstr(2:5);
        nstr{5} = ['There is not an equilibrium point near (',...
            num2str(z0(1)),', ', num2str(z0(2)), ').  Ready'];
        set(dud.notice,'string',nstr);
        return
    end
    zero =find(abs(z) < epsilon);
    z(zero) = zeros(size(zero));
    D=det(A); T=trace(A);
    if (D<-epsilon)
        string1 = ['There is a saddle point at (',...
            num2str(z(1)),', ', num2str(z(2)), ').'];
        string2 = '';
        type = 1;
    elseif(D>epsilon)
        if (T*T>4*D+epsilon)
            if (T<0)
                string1 = ['There is a nodal sink at (',...
                    num2str(z(1)),', ', num2str(z(2)), ').'];
                string2 = '';
                type = 2;
            else
                string1 = ['There is a nodal source at (',...
                    num2str(z(1)),', ', num2str(z(2)), ').'];
                string2 = '';
                type = 3;
            end
        elseif(T*T<4*D-epsilon)
            if(T<-epsilon)
                string1 = ['There is a spiral sink at (',...
                    num2str(z(1)),', ', num2str(z(2)), ').'];
                string2 = '';
                type = 4;
            elseif(T>epsilon)
                string1 = ['There is a spiral source at (',...
                    num2str(z(1)),', ', num2str(z(2)), ').'];
                string2 = '';
                type = 5;
            else
                string1 = ['There is a spiral equilibrium point at (',...
                    num2str(z(1)),', ', num2str(z(2)), ').'];
                string2 = 'Its specific type has not been determined.';
                type = 6;
            end
        else
            if (T>epsilon)
                string1 = ['There is a source at (',...
                    num2str(z(1)),', ', num2str(z(2)), ').'];
                string2 = 'Its specific type has not been determined.';
                type = 7;
            elseif (T<-epsilon)
                string1 = ['There is a sink at (',...
                    num2str(z(1)),', ', num2str(z(2)), ').'];
                string2 = 'Its specific type has not been determined.';
                type = 8;
            else
                string1 = ['There is an equilibrium point at (',...
                    num2str(z(1)),', ', num2str(z(2)), ').'];
                string2 = 'Its specific type has not been determined.';
                type = 9;
            end
        end
    else
        string1 = ['There is an equilibrium point at (',...
            num2str(z(1)),', ', num2str(z(2)), ')'];
        string2 = 'Its specific type has not been determined.';
        type = 9;
    end
    
    [V,D] = eig(A);
    EqPtList = dud.eqpts;
    infostr{1,1} = string1;
    infostr{2,1} = string2;
    infostr{3,1} = ' ';
    infostr{4,1} = 'The Jacobian is:';
    infostr{5,1} = ['      ',num2str(A(1,1)),'   ',num2str(A(1,2))];
    infostr{6,1} = ['      ',num2str(A(2,1)),'   ',num2str(A(2,2))];
    infostr{7,1} = '';
    infostr{8,1} = 'The eigenvalues and eigenvectors are:';
    infostr{9,1} = ['     ',num2str(D(1,1)),'    (',num2str(V(1,1)),', ',num2str(V(2,1)),')'];
    infostr{10,1} = ['     ',num2str(D(2,2)),'    (',num2str(V(1,2)),', ',num2str(V(2,2)),')'];
    
    k=1; l=size(EqPtList,1);
    while ((k <= l))
        if (norm((EqPtList(k,1:2)-z')./DY') <= 1e-3)
            break;
        end
        k = k+1;
    end
    if (k > l)
        EqPtList = [EqPtList;z',type];
        dud.eqpts = EqPtList;
        newh = plot(z(1),z(2),...
            'color',col,...
            'markersize',20,...
            'marker','.','Erasemode','none');
        dud.ephand = [dud.ephand;newh];
        set(newh,'Erasemode','normal');
        drawnow
    end
    ppeqpt = findobj('name','matpplane Equilibrium point data');
    
    if (isempty(ppeqpt))
        ppeqpt = figure('vis','off','resize','on',...
            'name','matpplane Equilibrium point data',...
            'numb','off','tag','matpplane');
        matpplane('figdefault',ppeqpt);
        set(ppeqpt,'menubar','none');
        
        ud.frame = uicontrol('style','frame');
        ud.eptext = uicontrol('style','text','string',infostr,'hor','left');
        ud.goaway = uicontrol('style','push','string','Go away',...
            'CallBack','close');
        ud.display = uicontrol('style','push',...
            'string','Display the linearization',...
            'CallBack','matpplane(''linear'')');
    else
        figure(ppeqpt);
        ud = get(ppeqpt,'UserData');
        set(ppeqpt,'vis','off');
        set(findobj(ppeqpt,'type','uicontrol'),'units','pixels');
        set(ud.eptext,'string',infostr);
    end
    ud.jac = A;
    ud.vectors = V;
    ud.type = type;
    ud.system = dud.syst;
    ud.settings = dud.settings;
    ud.color = dud.color;
    set(ppeqpt,'UserData',ud);
    left = 5; nudge = 3;
    ext = get(ud.eptext,'extent');
    n = size(infostr,1);
    ht = ext(4)/n+2*nudge;
    frbot = 3*left + 2*ht;
    txtbot = frbot + nudge;
    txtl = left + nudge;
    twind = [txtl,txtbot,ext(3)+nudge,ext(4)];
    frw = ext(3) + 3*nudge;
    frh = ext(4) + 2*nudge;
    frwind = [left,frbot,frw,frh];
    figw = frw + 2*left;
    figh = frh + 4*left + 2*ht;
    uni = get(0,'units');
    set(0,'units','pixels');
    ss = get(0,'screensize');
    set(0,'units',uni);
    sh = ss(4);
    figbot = sh - figh - 350;
    figwind = [30,figbot,figw,figh];
    buttw = frw - 4;
    buttl = left + 2;
    closewind = [buttl,left,buttw,ht];
    dispb = 2*left + ht;
    dispwind = [buttl,dispb,buttw,ht];
    set(ppeqpt,'pos',figwind);
    set(ud.frame,'pos',frwind);
    set(ud.eptext,'pos',twind);
    set(ud.goaway,'pos',closewind);
    set(ud.display,'pos',dispwind);
    set(findobj(ppeqpt,'type','uicontrol'),'units','normal');
    set(ppeqpt,'vis','on');
    set(get(ppeqpt,'child'),'vis','on');
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'Ready.';
    set(dud.notice,'string',nstr);
    set(ppdisp,'UserData',dud);
    
elseif strcmp(action,'stunst')
    
    ppdisp = findobj('name','matpplane Display');
    dud = get(ppdisp,'UserData');
    ecol = dud.color.eqpt;
    scol = dud.color.sep;
    ppdispa = dud.axes;
    aud = get(ppdispa,'UserData');
    DY = aud.DY;
    settings = dud.settings;
    dfcn = dud.function;
    EqPtList = dud.eqpts;
    Stop = norm(DY)*1e-4;
    
    %  Plot the stable and unstable orbits at a saddle point.
    
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'Choose a saddle point with the mouse.';
    set(dud.notice,'string',nstr);
    %    z0 = ppginput(1);
    z0 = ginput(1);
    %z = zeros(2,1); % unused preallocated variable
    z = matpplane('newton',z0,dfcn);
    
    flag = z(:,4);
    A = z(:,2:3);
    k = find(abs(A)<1e-8);
    A(k) = zeros(size(k));
    z = z(:,1);
    
    if ~flag | norm((z-z0')./DY)> 1/5 %#ok<OR2> % flag not scalar. can't use ||
        nstr(1:4) = nstr(2:5);
        nstr{5} = ['There is not an equilibrium point near (',...
            num2str(z0(1),2),', ', num2str(z0(2),2), ').'];
        set(dud.notice,'string',nstr);
        return
    end
    zero = find(abs(z) < 1e-4);
    z(zero) = zeros(size(zero));
    D=det(A);
    if (D>=0)
        nstr(1:4) = nstr(2:5);
        nstr{5} = ['The equilibrium point at (',...
            num2str(z(1),2),', ', num2str(z(2),2),...
            ') is not a saddle point.'];
        set(dud.notice,'string',nstr);
        return
    end
    if isempty(EqPtList)
        EqPtList = [z',1];
        dud.eqpts = EqPtList;
    else
        k=1; l=size(EqPtList,1);
        while ((k <= l))
            if (norm(EqPtList(k,1:2)-z') <= Stop)
                break;
            end
            k = k+1;
        end
        if (k > l)
            EqPtList = [EqPtList;z',1];
            dud.eqpts = EqPtList;
        end
    end
    nstr(1:4) = nstr(2:5);
    nstr{5} = ['Plotting from the saddle point at (',...
        num2str(z(1),2),', ', num2str(z(2),2),').'];
    set(dud.notice,'string',nstr);
    [V,L] = eig(A);
    [L,I] = sort(diag(L));
    
    %magn = settings.magn;
    refine = settings.refine;
    %ssize = settings.stepsize;
    %solver = settings.solver;
    tol = settings.tol;
    offset = norm(DY)/1000;
    lm=abs(L);
    lm=lm/max(lm);
    %lm = max([lm,0.2*ones(size(lm))]'))'; % Just wants the max in the 2nd dimension
    lm = max([lm,0.2*ones(size(lm))],[],2);
    offs = offset./lm;
    ud = get(dud.axes,'UserData');
    atol = tol*ud.DY*1e-4';
    %options = odeset('OutputFcn',@matppout,'Refine',refine,...
    %    'RelTol',tol,'Abstol',atol);
    stopbutt = findobj('tag','stop');
    set(stopbutt,'vis','on','enable','on');
    VV=V(:,I(1));
    newhand = zeros(4,1);
    w = z + offs(1)*VV;
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'The first stable trajectory';
    set(dud.notice,'string',nstr);
    int = [0,-1e6];
    
    solver = settings.solver;
    switch solver
        case 'Dormand Prince'
            solh = @matppdp45;
            opt = ppdisp;
        case 'Runge-Kutta 4'
            solh = @matpprk4;
            opt = ppdisp;
        case 'ode45'
            solh = @ode45;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',tol,'Abstol',atol);
        case 'ode23'
            solh = @ode23;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',tol,'Abstol',atol);
        case 'ode113'
            solh = @ode113;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',tol,'Abstol',atol);
        case 'ode15s'
            solh = @ode15s;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',tol,'Abstol',atol);
        case 'ode23s'
            solh = @ode23s;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',tol,'Abstol',atol);
        case 'ode23t'
            solh = @ode23t;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',tol,'Abstol',atol);
        case 'ode23tb'
            solh = @ode23tb;
            opt = odeset('OutputFcn',@matppout,'Refine',refine,...
                'RelTol',tol,'Abstol',atol);
    end
    %exist(dfcn);  % Does nothing
    dfh = str2func(dfcn);
    %cflag = 0;
    %[tp,Xp] = feval(solh,dfh,int,w,opt);
    [tp,Xp] = solh(dfh,int,w,opt);
    
    
    %   [tp,Xp] = ppdp45(dfcn,[0,-1e6],w,ppdisp);
    set(stopbutt,'enable','off');
    aud = get(ppdispa,'UserData');
    newhand(1) = aud.line;
    X1=[z';Xp];
    x1 = [X1,[NaN;tp]];
    set(newhand(1),'xdata',x1(:,1),'ydata',x1(:,2),'zdata',x1(:,3));
    
    w = z - offs(1)*VV;
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'The second stable trajectory';
    set(dud.notice,'string',nstr);
    set(stopbutt,'enable','on');
    %    [tp,Xp] = feval(solh,dfh,int,w,opt);
    [tp,Xp] = solh(dfh,int,w,opt);
    
    %   [tp,Xp] = ppdp45(dfcn,[0,-1e6],w,ppdisp);
    set(stopbutt,'enable','off');
    
    aud = get(ppdispa,'UserData');
    newhand(2) = aud.line;
    
    X2=[z';Xp];
    x2 = [X2,[NaN;tp]];
    set(newhand(2),'xdata',x2(:,1),'ydata',x2(:,2),'zdata',x2(:,3));
    
    VV=V(:,I(2));
    w = z + offs(2)*VV;
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'The first unstable trajectory';
    set(dud.notice,'string',nstr);
    set(stopbutt,'enable','on');
    int = [0,1e6];
    [tp,Xp] = solh(dfh,int,w,opt);
    set(stopbutt,'enable','off');
    
    aud = get(ppdispa,'UserData');
    newhand(3) = aud.line;
    X3=[z';Xp];
    x3 = [X3,[NaN;tp]];
    set(newhand(3),'xdata',x3(:,1),'ydata',x3(:,2),'zdata',x3(:,3));
    
    w = z - offs(2)*VV;
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'The second unstable trajectory';
    set(dud.notice,'string',nstr);
    set(stopbutt,'enable','on');
    [tp,Xp] = solh(dfh,int,w,opt);
    
    set(stopbutt,'vis','off');
    
    aud = get(ppdispa,'UserData');
    newhand(4) = aud.line;
    X4=[z';Xp];
    x4 = [X4,[NaN;tp]];
    set(newhand(4),'xdata',x4(:,1),'ydata',x4(:,2),'zdata',x4(:,3));
    
    eqpt = plot(z(1),z(2),...
        'color',ecol,...
        'markersize',20,...
        'marker','.','Erasemode','normal');
    dud.solhand = [dud.solhand;newhand(:)];
    dud.ephand = [dud.ephand;eqpt];
    set(newhand,'color',scol);
    set(newhand,'erase','normal');
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'Ready.';
    set(dud.notice,'string',nstr);
    set(gcf,'UserData',dud);
    
elseif strcmp(action,'zoomin')
    
    % 'zoomin' is the CallBackback for the Zoomin menu item.  It allows the
    % UserData to pick a new display rectangle.
    
    set(gcf,'WindowButtonDownFcn','matpplane(''zoom'')',...
        'WindowButtonUpFcn','1;','Interrupt','on');
    set(gca,'Interrupt','on');
    dud = get(gcf,'UserData');
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = ['Pick a new display rectangle by clicking and ',...
        'dragging the mouse, or by clicking on a point.'];
    set(dud.notice,'string',nstr);
    
elseif strcmp(action,'zoom')
    
    disph = gcf;
    dud = get(disph,'UserData');
    axh = dud.axes;
    aud = get(axh,'UserData');
    DY = aud.DY;
    w = dud.syst.wind;
    q1 = get(disph,'currentpoint');
    p1 = get(axh,'currentpoint');
    p1 = p1(1,1:2);
    rbbox([q1 0 0],q1);
    p2 = get(axh,'currentpoint');
    p2 = p2(1,1:2);
    if all(abs(p1'-p2')>0.01*DY)
        a = [p1;p2];
        a = [min(a);max(a)];
        DY = (a(2,:) - a(1,:))';
    else
        DY = DY/4;
        a(1) = max(w(1),p1(1)-DY(1));
        a(2) = min(w(2),p1(1)+DY(1));
        a(3) = max(w(3),p1(2)-DY(2));
        a(4) = min(w(4),p1(2)+DY(2));
        DY(1) = a(2) - a(1);
        DY(2) = a(4) - a(3);
    end
    WINvect = a(:)';
    dud.syst.wind = WINvect;
    aud.DY = DY;
    dwind = [WINvect(1); WINvect(3); -WINvect(2); -WINvect(4)];
    aud.cwind = dwind - dud.settings.magn*[aud.DY;aud.DY];
    set(axh,'UserData',aud);
    set(disph,'UserData',dud);
    set(disph,'WindowButtonDownFcn','matpplane(''down'')',...
        'WindowButtonUpFcn','');
    matpplane('dirfield',disph);
    ppset = findobj('name','matpplane Setup');
    if isempty(ppset)
        matpplane('confused');
    else
        sud = get(ppset,'UserData');
        sud.c.wind = WINvect;
        sud.o.wind = WINvect;
        set(sud.h.wind(1),'string',num2str(WINvect(1)));
        set(sud.h.wind(2),'string',num2str(WINvect(2)));
        set(sud.h.wind(3),'string',num2str(WINvect(3)));
        set(sud.h.wind(4),'string',num2str(WINvect(4)));
        set(ppset,'UserData',sud);
    end
    
elseif strcmp(action,'showbar')
    
    sbfig = gcbf;
    domymenu('menubar','toggletoolbar',sbfig);
    hhset = get(0,'showhiddenhandles');
    set(0,'showhiddenhandles','on');
    state = get(sbfig,'toolbar');
    if strcmp(state,'figure')
        name = get(sbfig,'name');
        fixtb = 'set(gcbo,''state'',''off'');';
        set(findobj(sbfig,'tooltipstr','Print'),...
            'clickedCallBackback','matpplane(''print'');');
        switch name
            case 'matpplane Display'
                set(findobj(sbfig,'tooltipstr','Zoom Out'),...
                    'clickedCallBackback',['matpplane(''zoomback'');' fixtb]);
                set(findobj(sbfig,'tooltipstr','Zoom In'),...
                    'clickedCallBackback',['matpplane(''zoomin'');' fixtb]);
                delete(findobj(sbfig,'tooltipstr','Rotate 3D'));
                
            case 'matpplane Linearization'
                delete(findobj(sbfig,'tooltipstr','Zoom Out'));
                delete(findobj(sbfig,'tooltipstr','Zoom In'));
                delete(findobj(sbfig,'tooltipstr','Rotate 3D'));
            otherwise
                delete(findobj(sbfig,'tooltipstr','Zoom Out'));
                delete(findobj(sbfig,'tooltipstr','Zoom In'));
                fud = get(sbfig,'UserData');
                switch fud.type
                    case {4, 5}
                        set(findobj(sbfig,'tooltipstr','Rotate 3D'),'vis','on');
                    otherwise
                        set(findobj(sbfig,'tooltipstr','Rotate 3D'),'vis','off');
                end
                
        end
        sbmh = findobj(sbfig,'label','Show &Toolbar');
        set(sbmh,'label','Hide &Toolbar','checked','off');
        uni = get(0,'units');
        ss = get(0,'screensize');
        funit = get(sbfig,'units');
        set(sbfig,'units',uni);
        %sw = ss(3);
        sh = ss(4);
        fpos = get(sbfig,'pos');
        if fpos(2)+fpos(4)>sh-40
            fpos(2) = sh - fpos(4) -70;
            set(sbfig,'pos',fpos);
        end
        set(sbfig,'units',funit);
        
    else
        sbmh = findobj(sbfig,'label','Hide &Toolbar');
        set(sbmh,'label','Show &Toolbar','checked','off');
    end
    set(0,'showhiddenhandles',hhset)
    
elseif strcmp(action,'dall')
    
    % 'dall' is the CallBackback for the Erase all graphics objects.
    
    disph = gcf;
    dud = get(disph,'UserData');
    notice = dud.notice;
    kk = find(ishandle(dud.solhand));
    delete(dud.solhand(kk));
    dud.solhand = [];
    kk = find(ishandle(dud.ephand));
    delete(dud.ephand(kk));
    dud.ephand = [];
    dud.eqpts = [];
    kk = find(ishandle(dud.contours));
    delete(dud.contours(kk));
    dud.contours = [];
    kk = find(ishandle(dud.ics));
    delete(dud.ics(kk));
    dud.ics = [];
    if ishghandle(notice)
        set(dud.butt(1),'enable','off');
    end
    set(disph,'UserData',dud);
    
    
elseif strcmp(action,'dallsol')
    
    % 'dallsol' is the CallBackback for the Erase all solutions option.
    
    disph = gcf;
    dud = get(disph,'UserData');
    notice = dud.notice;
    kk = find(ishandle(dud.solhand));
    delete(dud.solhand(kk));
    dud.solhand = [];
    if ishghandle(notice)
        set(dud.butt(1),'enable','off');
    end
    set(disph,'UserData',dud);
    
elseif strcmp(action,'dallep')
    
    % 'dallep' is the CallBackback for the Erase all equilibrium points option.
    
    disph = gcf;
    dud = get(disph,'UserData');
    notice = dud.notice;
    kk = find(ishandle(dud.ephand));
    delete(dud.ephand(kk));
    dud.ephand = [];
    dud.eqpts = [];
    if ishghandle(notice)
        set(dud.butt(1),'enable','off');
    end
    set(disph,'UserData',dud);
    
elseif strcmp(action,'dalllev')
    
    % 'dalllev' is the CallBackback for the Erase all level curves option.
    
    disph = gcf;
    dud = get(disph,'UserData');
    notice = dud.notice;
    kk = find(ishandle(dud.contours));
    delete(dud.contours(kk));
    dud.contours = [];
    if ishghandle(notice)
        set(dud.butt(1),'enable','off');
    end
    set(disph,'UserData',dud);
    
elseif strcmp(action,'dallics')
    
    % 'dallics' is the CallBackback for the Erase all marked initial cond.
    
    disph = gcf;
    dud = get(disph,'UserData');
    notice = dud.notice;
    kk = find(ishandle(dud.ics));
    delete(dud.ics(kk));
    dud.ics = [];
    if ishghandle(notice)
        set(dud.butt(1),'enable','off');
    end
    set(disph,'UserData',dud);
    
elseif strcmp(action,'newton')
    
    % Newton's method to find equilibrium points.
    
    Iterlimit = 50;
    
    Iter=0;
    
    % The increment for calculating approximate derivatives.
    
    h=.000001;
    
    dfcn = input2;
    
    zInit = input1;
    
    zNext=zInit(:);
    % functionf = zeros(2,1); % Unused preallocated variable
    functionf=feval(dfcn,0,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.
        Jacobian=zeros(2,2);
        for jjw=1:2
            sav = zNext(jjw);
            zNext(jjw) = zNext(jjw) + h;
            functionfhj = feval(dfcn,0,zNext);
            functionfhj = functionfhj(:);
            Jacobian(:,jjw) = (functionfhj-functionf)/h;
            zNext(jjw) = sav;
        end
        zNext = zNext - Jacobian\(functionf);
        functionf = feval(dfcn,0,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(dfcn,0,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 = deblank(input1);
    value = input2;
    value = ['(',value,')'];
    str = input3;
    ll = length(str);
    lp = length(para);
    %lv = length(value);
    
    if strcmp(para,str)
        str = value;
    elseif (ll >= lp+1)
        k = strfind(str,para);
        
        lk = length(k);
        lopstr = '(+-*/^';
        ropstr = ')+-*/^';
        s = [];
        pos = 1;
        for jj = 1:lk
            if ((k(jj) == 1)||(find(lopstr == str(k(jj)-1))))...
                    &&((k(jj)+lp-1 == ll)||(find(ropstr == str(k(jj) + lp))))
                s = [s,str(pos:(k(jj)-1)),value]; %#ok<AGROW>
                pos = k(jj)+lp;
            end
        end
        str = [s,str(pos:ll)];
    end
    output = str;
    
    
elseif strcmp(action,'settings')
    
    % 'settings' is the CallBack back for the Settings menu option.  It sets
    % up the matpplane Settings window, which allows the UserData to interactively
    % change several parameters that govern the behaviour of the program.
    
    % This gets called when user clicks "Settings" or "Solver" in
    % Display->Options
    
    dud = get(gcf,'UserData');
    data.settings = dud.settings;
    solver = dud.settings.solver;
    left = 2; nudge = 3; varl = 60;
    setfig = findobj('name','matpplane Settings');
    if ~isempty(setfig)
        delete(setfig);
    end
    
    setfig = figure('name','matpplane Settings',...
        'numb','off',...
        'tag','matpplane','vis','off');
    
    matpplane('figdefault',setfig);
    set(setfig,'menubar','none');
    
    setfr = uicontrol('style','frame');
    speedfr = uicontrol('style','frame');
    
    cwfr = uicontrol('style','frame');
    ss = uicontrol('style','text','horiz','center',...
        'string',['Settings for ',dud.settings.solver]);
    
    nstepCallBack =[
        'data = get(gcf,''UserData'');',...
        'ss = max(round(str2num(get(data.nstep,''string''))),0);',...
        'data.settings.refine = ss;',...
        'set(data.nstep,''string'',num2str(ss));',...
        'set(gcf,''UserData'',data);'];
    
    nstep = uicontrol('style','text','horiz','left',...
        'string','Number of plot steps per computation step:  ');
    data.nstep = uicontrol('style','edit','CallBack',nstepCallBack,...
        'string',num2str(data.settings.refine));
    
    if strcmp(solver,'Runge-Kutta 4')
        rtolstr = 'Step size: ';
        rtolentry = num2str(data.settings.stepsize);
        rtolCallBack =['data = get(gcf,''UserData'');',...
            'data.settings.stepsize = str2num(get(data.rtol,''string''));',...
            'set(data.rtol,''string'',num2str(data.settings.stepsize));',...
            'set(gcf,''UserData'',data);'];
    else
        rtolstr = 'Relative error tolerance: ';
        rtolentry = num2str(data.settings.tol);
        rtolCallBack =['data = get(gcf,''UserData'');',...
            'data.settings.tol = str2num(get(data.rtol,''string''));',...
            'set(data.rtol,''string'',num2str(data.settings.tol));',...
            'set(gcf,''UserData'',data);'];
    end
    rtol = uicontrol('style','text','horiz','left',...
        'string',rtolstr);
    data.rtol = uicontrol('style','edit','CallBack',rtolCallBack,...
        'string',rtolentry);
    
    kk = 1+2*data.settings.magn;
    magCallBack = ['data = get(gcf,''UserData'');',...
        'mag = 	(str2num(get(data.mag,''string''))-1)/2;',...
        'data.settings.magn = max(0,mag);',...
        'set(gcf,''UserData'',data);'];
    
    speedCallBack = ['data = get(gcf,''UserData'');',...
        'me = data.speed;',...
        'val = round(get(me,''value''));',...
        'set(me,''value'',val);',...
        'set(data.sp.val,''string'',num2str(val));',...
        'data.settings.speed = val;',...
        'set(gcf,''UserData'',data);'];
    data.speed = uicontrol('style','slider',...
        'string','Steps per second.',...
        'min',2,...
        'max',100,...
        'value',data.settings.speed,...
        'sliderstep',[1/98,10/98],...
        'CallBack',speedCallBack);
    
    data.sp.min = uicontrol('style','text','string',' 2',...
        'horiz','left');
    data.sp.max = uicontrol('style','text','string','100 ',...
        'horiz','right');
    data.sp.val = uicontrol('style','text',...
        'string',num2str(data.settings.speed),...
        'horiz','center');
    
    pps = uicontrol('style','text',...
        'string','Solution steps per second.');
    
    
    
    
    cw1 = uicontrol('style','text','horiz','left',...
        'string','The calculation window is');
    data.mag = uicontrol('style','edit','CallBack',magCallBack,...
        'string',num2str(kk));
    cw2 = uicontrol('style','text','horiz','left',...
        'string',' times as large as the ');
    cw3 = uicontrol('style','text','horiz','left',...
        'string','display window.');
    
    gob = uicontrol('style','push','string','Go Away','CallBack','close');
    
    chb = uicontrol('style','push',...
        'string','Change settings',...
        'CallBack','matpplane(''setchange'')');
    
    frsep = 1;
    ext = get(nstep,'extent');
    ht = ext(4)+nudge;
    stw = ext(3);		% = nstepw = rtolw
    varl = varl*ht/19;
    cwfrb = 2*left + ht;
    cwfrh = 2*nudge + 2*ht;
    speedfrb = cwfrb + cwfrh + frsep;
    speedfrh = 4*nudge + 3*ht;
    setfrb = speedfrb + speedfrh + frsep;
    setfrh = 3*nudge + 3*ht;
    figh = setfrb + setfrh + left;
    bb = left;
    cw3b = cwfrb + nudge;
    cw2b = cw3b + ht;	% = cw1b = cweb
    rtolb = setfrb + nudge;
    rtolbb = rtolb;
    if strcmp(solver,'Runge-Kutta 4')
        rtolbb = rtolbb +ht/2;
    end
    nstepb = rtolb + ht;
    ssb = nstepb + ht + nudge;
    ssw = stw + varl;
    spb1 = speedfrb+2*nudge;
    spb2 = spb1+ht;
    spb3 = spb2+ht+nudge;
    sptw = varl;
    sl = left + nudge;
    sptl1 = sl;
    sptsep = (ssw-3*sptw-sl)/2;
    sptl2 = sptl1 + sptw + sptsep;
    sptl3 = sptl2 + sptw + sptsep;
    
    cw1ext = get(cw1,'extent');
    cw1w = cw1ext(3);
    cw2ext = get(cw2,'extent');
    cw2w = cw2ext(3);
    cwew = 40*ht/19;
    frw = 2*nudge + ssw;
    figw = frw + 2*left;
    buttw = figw/3;
    sep = figw/9;
    sl = left + nudge;
    gobl = sep;
    chbl = 2*sep + buttw;
    
    %sunit = get(0,'units');
    set(0,'units','pix');
    ssize = get(0,'screensize');
    figb = ssize(4) - figh - 50;
    
    set(setfig,'pos',[20,figb,figw,figh]);
    set(setfr,'pos',[left,setfrb,frw,setfrh]);
    set(speedfr,'pos',[left,speedfrb,frw,speedfrh]);
    set(cwfr,'pos',[left,cwfrb,frw,cwfrh]);
    set(ss,'pos',[sl,ssb,ssw,ht]);
    set(nstep,'pos',[sl,nstepb,stw,ht]);
    set(data.nstep,'pos',[sl+stw,nstepb,varl,ht]);
    set(rtol,'pos',[sl,rtolbb,stw,ht]);
    set(data.rtol,'pos',[sl+stw,rtolbb,varl,ht]);
    
    set(data.speed,'pos',[sl,spb1,ssw,ht],...
        'backgroundcolor',0.6*[1 1 1],...
        'foregroundcolor','r');
    set(data.sp.min,'pos',[sptl1,spb2,sptw,ht]);
    set(data.sp.max,'pos',[sptl3,spb2,sptw,ht]);
    set(data.sp.val,'pos',[sptl2,spb2,sptw,ht]);
    set(pps,'pos',[sl,spb3,ssw,ht]);
    
    set(cw1,'pos',[sl,cw2b,cw1w,ht]);
    set(cw2,'pos',[sl+cw1w+cwew,cw2b,cw2w,ht]);
    set(cw3,'pos',[sl,cw3b,ssw,ht]);
    set(data.mag,'pos',[sl+cw1w,cw2b,cwew,ht]);
    set(gob,'pos',[gobl,bb,buttw,ht]);
    set(chb,'pos',[chbl,bb,buttw,ht]);
    
    
    set(setfig,'UserData',data);
    set(setfig,'units','normal');
    set(findobj(setfig,'type','uicontrol'),'units','normal');
    set(setfig,'vis','on','resize','on');
    ch = get(setfig,'child');
    set([nstep,data.nstep],'vis','off')
    if strcmp(solver,'Runge-Kutta 4')
        ch = setdiff(ch,[nstep,data.nstep]);
    end
    edith = findobj(setfig,'style','edit');
    set(edith,'backgroundcolor','w');
    set(ch,'vis','on');
    
    
elseif strcmp(action,'setchange')
    
    % 'setchange' is the CallBackback for the Change button
    % on the matpplane Settings window.
    
    data = get(gcf,'UserData');
    settings = data.settings;
    ppdisp = findobj('name','matpplane Display');
    dud = get(ppdisp,'UserData');
    if isempty(ppdisp)
        matpplane('confused');
    else
        dud.settings = settings;
        set(ppdisp,'UserData',dud);
        WINvect = dud.wmat;
        WINvect = WINvect(size(WINvect,1),:);
        dwind = [WINvect(1); WINvect(3); -WINvect(2); -WINvect(4)];
        aud = get(dud.axes,'UserData');
        aud.cwind = dwind - dud.settings.magn*[aud.DY;aud.DY];
        set(dud.axes,'UserData',aud);
    end
    ppset = findobj('name','matpplane Setup');
    if isempty(ppset)
        matpplane('confused');
    else
        sud = get(ppset,'UserData');
        sud.settings = settings;
        set(ppset,'UserData',sud);
    end
    close
    
    
elseif strcmp(action,'delete')
    
    % 'delete' is the CallBackback for the Delete a graphics object selection
    % on the menu.
    
    disph = gcf;
    dud = get(disph,'UserData');
    arr = dud.arr;
    lv = get(arr.lines,'vis');
    av = get(arr.arrows,'vis');
    if ~isempty(arr.hx)
        nv = get(arr.hx(1),'vis');
    elseif ~isempty(arr.hy)
        nv = get(arr.hx(1),'vis');
    else
        nv = zeros(1,0);
    end
    if ~isempty(arr.barrows)
        bv = get(arr.barrows,'vis');
    else
        bv = zeros(1,0);
    end
    handles = [arr.lines;arr.arrows;arr.hx;arr.hy;arr.barrows];
    set(handles,'vis','off');
    %oldCallBack = get(disph,'WindowButtonDownFcn');
    set(disph,'WindowButtonDownFcn','');
    %trjh = dud.solhand;
    notice = dud.notice;
    if ishghandle(notice)   % Display window
        nstr = get(notice,'string');
        nstr(1:4) = nstr(2:5);
        nstr{5} = 'Select a graphics object with the mouse.';
        set(notice,'string',nstr);
    end
    %   ppginput(1);
    ginput(1);
    objh = get(disph,'currentobject');
    typ = get(objh,'type');
    axh = dud.axes;
    hh = get(dud.title.axes,'children');
    hh = [hh;get(axh,'title');get(axh,'xlabel');get(axh,'ylabel')];
    if ishghandle(notice)   % Display window
        eph = dud.ephand;
        %levh = dud.contours;
        eqpt = dud.eqpts;
        %ics = dud.ics;
        nstr(1:4) = nstr(2:5);
        if strcmp(typ,'line')
            eqk = find(eph == objh,1);
            if ~isempty(eqk)  % An equilibrium point.
                if size(eqpt,1) == 1
                    dud.eqpts = [];
                else
                    dud.eqpts = [eqpt(1:eqk-1,:);eqpt(eqk+1:size(eqpt,1),:)];
                end
            end
            dud.ics = setdiff(dud.ics,objh);
            dud.solhand = setdiff(dud.solhand,objh);
            dud.ephand = setdiff(dud.ephand,objh);
            dud.contours = setdiff(dud.contours,objh);
            delete(objh);
            nstr{5} = 'Ready.';
            
        elseif strcmp(typ,'text') && ~ismember(objh,hh)
            dud.contours = setdiff(dud.contours,objh)';
            delete(objh);
            nstr{5} = 'Ready.';
            
        else
            nstr{5} = 'The object you selected cannot be deleted.';
            
        end
        set(notice,'string',nstr);
    else     % Linearization window
        if strcmp(typ,'line')
            eqk = find(dud.ephand == objh,1);
            if isempty(eqk)  % Not an equilibrium point.
                dud.solhand = setdiff(dud.solhand,objh)';
                dud.ics = setdiff(dud.ics,objh);
                delete(objh);
            end
        elseif ~ismember(objh,hh)
            delete(objh);
        end
    end
    
    set(arr.lines,'vis',lv);
    set(arr.arrows,'vis',av);
    set([arr.hx;arr.hy],'vis',nv);
    set(arr.barrows,'vis',bv);
    set(disph,'UserData',dud);
    set(disph,'WindowButtonDownFcn','matpplane(''down'')');
    
elseif strcmp(action,'text')
    
    ppdisp = gcf;
    oldCallBack = get(ppdisp,'WindowButtonDownFcn');
    set(ppdisp,'WindowButtonDownFcn','');
    prompt = ['Enter the text here. Then choose ',...
        'the location in the Display Window.'];
    txtstr = inputdlg(prompt,'Text entry');
    if ~isempty(txtstr)
        txtstr = txtstr{1};
        figure(ppdisp);
        gtext(txtstr);
    end
    set(ppdisp,'WindowButtonDownFcn',oldCallBack);
    
elseif strcmp(action,'plotxy')  % Start a graph
    
    type = input1;
    
    ppdisp =gcf;
    dud = get(ppdisp,'UserData');
    
    handles = dud.solhand;
    %printstr = dud.printstr;
    exsol=0;
    lll = length(handles);
    
    if (lll)  % Are there any orbits?
        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)
        nstr = get(dud.notice,'string');
        nstr(1:3) = nstr(3:5);
        nstr{4} = 'There are no solution curves.';
        nstr{5} = 'Ready.';
        set(dud.notice,'string',nstr);
        return
    end
    
    oldCallBack = get(ppdisp,'WindowButtonDownFcn');
    set(ppdisp,'WindowButtonDownFcn','');
    arr = dud.arr;
    lv = get(arr.lines,'vis');
    av = get(arr.arrows,'vis');
    if ~isempty(arr.hx)
        nv = get(arr.hx(1),'vis');
    elseif ~isempty(arr.hy)
        nv = get(arr.hy(1),'vis');
    else
        nv = zeros(1,0);
    end
    if ~isempty(arr.barrows)
        bv = get(arr.barrows,'vis');
    else
        bv = zeros(1,0);
    end
    handles = [arr.lines;arr.arrows;arr.hx;arr.hy;arr.barrows];
    set(handles,'vis','off');
    
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'Select a solution curve with the mouse.';
    set(dud.notice,'string',nstr);
    % ppginput(1);
    ginput(1);
    objh = get(ppdisp,'currentobject');
    if isempty(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
        nstr = get(dud.notice,'string');
        nstr(1:4) = nstr(2:5);
        nstr{5} = 'This object is not a solution curve.';
        set(dud.notice,'string',nstr);
        set(ppdisp,'WindowButtonDownFcn',oldCallBack);
        set(arr.lines,'vis',lv);
        set(arr.arrows,'vis',av);
        set([arr.hx;arr.hy],'vis',nv);
        return
    end
    set(arr.lines,'vis',lv);
    set(arr.arrows,'vis',av);
    set([arr.hx;arr.hy],'vis',nv);
    set(arr.barrows,'vis',bv);
    t=get(objh,'zdata');
    x=get(objh,'xdata');
    y=get(objh,'ydata');
    hh = matpplane('plotxyfig',type,gcf);
    set(ppdisp,'WindowButtonDownFcn',oldCallBack);
    aud = get(hh(2),'UserData');
    ud = get(hh(1),'UserData');
    ud.data = [t;x;y];
    set(hh(1),'UserData',ud);
    matpplane('plxy',aud.rad,hh(1));
    
elseif strcmp(action,'plotxyfig')  % Build the graph figure.
    
    type = input1;
    cfig = input2;  % This is ppdisp or a ppxy
    ppdisp = findobj('name','matpplane Display');
    dud = get(ppdisp,'UserData');
    ppdispa = dud.axes;
    xstring = get(get(ppdispa,'XLabel'),'string');
    ystring = get(get(ppdispa,'YLabel'),'string');
    dud = get(ppdisp,'UserData');
    
    ppxy = figure('numbertitle','off',...
        'tag','matpplane','visible','off');
    name = ['matpplane t-plot # ',int2str(get(ppxy,'Number'))];
    ud.color = dud.color;
    ud.type = 0;
    set(ppxy,'name',name,'UserData',ud);
    matpplane('figdefault',ppxy);
    
    % Menus and Toolbar
    
    hhsetup = get(0,'showhiddenhandles');
    set(0,'showhiddenhandles','on');
    
    % Configure the Toolbar.
    
    set(ppxy,'ToolBar','none');
    
    % Menus
    
    %emenu = findobj(gcf,'label','&Edit');
    tmenu = findobj(gcf,'label','&Tools');
    delete(tmenu)
    
    % File menu
    
    fmenu = findobj(gcf,'label','&File');
    delete(findobj(fmenu,'label','&New Figure'));
    delete(findobj(fmenu,'label','&Open...'));
    set(findobj(fmenu,'label','&Close'),'pos',1);
    set(findobj(fmenu,'label','&Save'),...
        'pos',2,'separator','off');
    set(findobj(fmenu,'label','Save &As...'),...
        'pos',3);
    set(findobj(fmenu,'label','&Export...'),...
        'pos',4);
    delete(findobj(fmenu,'label','Pre&ferences...'));
    set(findobj(fmenu,'label','Pa&ge Setup...'),'pos',4);
    set(findobj(fmenu,'label','Print Set&up...'),'pos',5);
    set(findobj(fmenu,'label','Print Pre&view...'),'pos',6);
    set(findobj(fmenu,'label','&Print...'),'pos',7);
    [~] = uimenu(fmenu,'label','Quit matpplane',...
        'CallBack','matpplane(''quit'')',...
        'separator','on');
    
    % Insert menu
    
    imenu = findobj(gcf,'label','&Insert');
    inschild = get(imenu,'child');
    legitem = findobj(inschild,'label','&Legend');
    colitem = findobj(inschild,'label','&Colorbar');
    delete([legitem,colitem]);
    
    % View menu
    
    set(findobj(gcf,'label','&Figure Toolbar'),...
        'CallBack','matpplane(''showbar'')');
    
    set(0,'showhiddenhandles',hhsetup);
    
    fs = get(ppxy,'defaultaxesfontsize');
    r = fs/10;
    axw = 437*.8;	% Default axes width
    axh = 315*.8;	% Default axes height
    axl = 45;	% Default axes left
    %buttw = 40;     % Default button width
    titleh = 45;    % Default title height
    legw = 90;     % Default legend width
    axb = 35;       % Default axes bottom
    sep = 10;
    fh = axb + axh + titleh;  % Default figure height.
    fw = axl + axw + sep + legw + sep; % Default figure width.
    uni = get(0,'units');
    set(0,'units','pixels');
    ss = get(0,'screensize');
    set(0,'units',uni);
    sw = ss(3);sh = ss(4);
    if r*fh > sh -80
        r = (sh-80)/fh;
    end
    if r*fw > sw - 50
        r = (sw - 50)/fw;
    end
    fs = 10*r;
    lw = 0.5*r;
    set(ppxy,'defaultaxesfontsize',fs,'defaultaxeslinewidth',lw);
    set(ppxy,'defaulttextfontsize',fs);
    set(ppxy,'defaultlinelinewidth',lw);
    set(ppxy,'defaultuicontrolfontsize',fs*0.9);
    axw = r*axw;
    axh = r*axh;
    axl = r*axl;
    legw = r*legw;
    axb = r*axb;
    sep = r*sep;
    
    % The legend.
    
    
    leg = axes('units','pix',...
        'box','on',...
        'drawmode','fast',...
        'nextplot','add',...
        'xtick',-1,'ytick',-1,...
        'xticklabel','','yticklabel','',...
        'xlim',[0,1],'ylim',[0,1],...
        'clipping','on','visible','off');
    
    set(leg,'UserData',char(xstring,ystring));
    axes(leg);
    xh = text(0,0,xstring,'units','norm','visible','off',...
        'parent',leg);
    yh = text(0,0,ystring,'units','norm','visible','off',...
        'parent',leg);
    set(xh,'units','pix')
    set(yh,'units','pix')
    xext = get(xh,'extent');
    yext = get(yh,'extent');
    set(xh,'units','norm')
    set(yh,'units','norm')
    th = max(xext(4),yext(4))+3*r;
    axhn = max(axh, 11*th + 4 + 6*r);  % Nudge the axes height if needed.
    axw = axw*axhn/axh;
    axh = axhn;
    frsep = (axh - 11*th - 4 - 6*r)/4;
    
    legh = 2*(th+3*r);
    legb = axb + axh - legh;  %frsep*4 + th*9 +4;
    legl = axl + axw + sep;
    
    
    % If the legend text is too big, make the legend and the figure larger.
    
    legpos = [legl+1, legb, legw, legh];
    set(leg,'pos',legpos);
    xext = get(xh,'extent');
    yext = get(yh,'extent');
    tw = max(xext(3),yext(3))+0.1;
    
    rrrr = 0.5;
    if (tw > rrrr)
        rrr = tw/rrrr;
        legw = legw*rrr;
        tw = rrrr;
    end
    legpos = [legl+1, legb, legw, legh];
    set(leg,'pos',legpos);
    fw = axl + axw + sep + legw + sep; % Figure width.
    
    tx = min(1-tw,2/3);
    set(xh,'pos',[tx,2/3]);
    set(yh,'pos',[tx,1/3]);
    
    line('xdata',[0.1,tx-0.1],'ydata',[2/3,2/3],'linestyle','-',...
        'color',ud.color.tx,'vis','off');
    line('xdata',[0.1,tx-0.1],'ydata',[1/3,1/3],'linestyle','--',...
        'color',ud.color.ty,'vis','off');
    
    % Set up the title.
    
    tax = axes('box','off',...
        'xlim',[0 1],'ylim',[0 1],...
        'units','pix','vis','off',...
        'xtick',-1,'ytick',-1,...
        'xticklabel','','yticklabel','');
    titb = axb + axh;
    titl = axl;
    titw = axw;
    titeq = text(0.01,0.5,dud.tstr,'vert','middle');
    set(titeq,'units','pix');
    ext = get(titeq ,'extent');
    set(titeq,'units','data');
    titleh = ext(4)+6*r;
    p1 = get(dud.title.p1,'string');
    p2 = get(dud.title.p2,'string');
    p3 = get(dud.title.p3,'string');
    pp1 = text(0.85,0.5,p1,'vert','middle');
    ext = get(pp1,'extent');
    pp1l = 1 - ext(3);
    set(pp1,'pos',[pp1l,0.5]);
    pp2 = text(0.7,0.5,p2,'vert','middle');
    ext = get(pp2,'extent');
    pp2l = pp1l - ext(3)-0.05;
    set(pp2,'pos',[pp2l,0.5]);
    pp3 = text(0.7,0.5,p3,'vert','middle');
    ext = get(pp3,'extent');
    pp3l = pp2l - ext(3)-0.05;
    set(pp3,'pos',[pp3l,0.5]);
    
    fh = axb + axh + titleh;  % Figure height.
    
    pos = get(cfig,'pos');
    fl = max(0,pos(1)-30);
    fb = max(20,pos(2)-30);
    fpos = [fl,fb,fw,fh];
    set(ppxy,'pos',fpos);
    
    set(tax,'pos',[titl,titb,titw,titleh]);
    
    set(tax,'color',get(ppxy,'color'));
    
    
    % The t-plot axes.
    
    axpos = [axl,axb,axw,axh];
    
    axy = axes('units','pix','pos',axpos,'box',...
        'on','xgrid','on','ygrid','on','next','add',...
        'drawmode','fast');
    
    output = [ppxy,axy];
    
    
    
    % The position of the close button
    
    bbot = axb;
    bwid = legw;
    bleft = legl;
    bpos = [bleft,bbot,bwid,th];
    
    % The position of the print button.
    
    pbbot = bbot + th + frsep;
    pbpos = [bleft, pbbot, legw, th];
    
    % The position of the crop button.
    cbbot = pbbot + th + frsep;
    cbpos = [bleft, cbbot, legw, th];
    
    % The positions of the radio frame and its elements.
    
    frbot = cbbot + th + frsep;
    frleft = legl;
    frht = 6*th+4;
    frpos = [frleft,frbot,legw,frht];
    
    frtitleft = frleft+2;
    frtitwid = legw-4;
    frtitpos = [frtitleft, frbot+5*th, frtitwid, th];
    
    radleft = frleft + 2;
    radwid = legw - 4;
    radbot5 = frbot + 2;
    radbot4 = radbot5 + th;
    radbot3 = radbot4 + th;
    radbot2 = radbot3 + th;
    radbot1 = radbot2 + th;
    pCallBack = [
        'ppdisp = findobj(''name'',''matpplane Display'');',...
        'dud = get(ppdisp,''UserData'');',...
        'eval(dud.printstr)'
        ];
    
    [~] = uicontrol('style','push',...
        'pos',bpos,...
        'string','Go away',...
        'CallBack','close');
    
    [~] = uicontrol('style','push',...
        'pos',pbpos,...
        'string','Print',...
        'CallBack',pCallBack);
    
    cbut = uicontrol('style','push',...
        'pos',cbpos,...
        'string','Crop',...
        'CallBack','matpplane(''crop'')');
    
    [~] = uicontrol('style','frame','pos',frpos);
    
    [~] = 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 th],...
        'string',xstring,'value',rval1);
    
    rad(2) = uicontrol('style','radio',...
        'pos',[radleft radbot2 radwid th],...
        'string',ystring,'value',rval2,'max',2);
    
    rad(3) = uicontrol('style','radio',...
        'pos',[radleft radbot3 radwid th],...
        'string','Both','value',rval3,'max',3);
    
    rad(4) = uicontrol('style','radio',...
        'pos',[radleft radbot4 radwid th],...
        'string','3 D','value',rval4,'max',4);
    
    rad(5) = uicontrol('style','radio',...
        'pos',[radleft radbot5 radwid th],...
        'string','Composite','value',rval5,'max',5);
    
    for i=1:5
        set(rad(i),'UserData',[rad(:,[1:(i-1),(i+1):5]),leg,axy]);
    end
    
    CallBackrad = [
        'me = gcbo;',...
        'vv = get(me,''value'');'...
        'mm = get(me,''max'');'...
        'if vv ,'...
        '	hand = get(me,''UserData'');',...
        '	set(hand(1:4),''value'',0);',...
        '	matpplane(''plxy'',me,gcf);'...
        'end, '...
        'set(me,''value'',mm);'...
        'axy = gca;',...
        'aud = get(axy,''UserData'');',...
        'aud.rad = me;',...
        'set(axy,''UserData'',aud);'];
    
    
    set(rad,'CallBack',CallBackrad);
    
    set(findobj(ppxy,'type','axes'),'units','normal');
    set(findobj(ppxy,'type','uicontrol'),'units','normal');
    set(ppxy,'resize','on');
    
    set(ppxy,'visible','on');
    
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'Ready.';
    set(dud.notice,'string',nstr);
    aud.h = [];
    aud.int = [0 0];
    aud.crop = cbut;
    aud.rad = rad(type);
    aud.leg = leg;
    set(axy,'UserData',aud);
    set(cbut,'UserData',axy);
    
elseif strcmp(action,'plxy')
    
    radbut = input1;
    fig = input2;
    fud = get(fig,'UserData');
    figure(fig)
    tbh = findobj(allchild(fig),'flat','type','uitoolbar');
    r3dh = findobj(tbh,'tooltipstr','Rotate 3D');
    ud = get(radbut,'UserData');
    axy = ud(6);
    axis('auto');
    delete(get(axy,'children'));
    leg = ud(5);
    legch = get(leg,'children');
    type = get(radbut,'max');
    fud.type = type;
    set(fig,'UserData',fud);
    if type == 3
        set(leg,'visible','on');
        set(legch,'visible','on');
    else
        set(leg,'visible','off');
        set(legch,'visible','off');
    end
    
    aud = get(axy,'UserData');
    aud.h = [];
    set(axy,'UserData',aud);
    strings = get(leg,'UserData');
    xstring = deblank(strings(1,:));
    ystring = deblank(strings(2,:));
    ud = get(fig,'UserData');
    data = ud.data;
    color = ud.color;
    
    t = data(1,:);
    x = data(2,:);
    y = data(3,:);
    tmin = min(t);
    tmax = max(t);
    xmin = min(x);
    xmax = max(x);
    ymin = min(y);
    ymax = max(y);
    et = max((tmax-tmin)/20,1e-4);
    ex = max((xmax-xmin)/20,1e-4);
    ey = max((ymax-ymin)/20,1e-4);
    tmin1 = tmin - et;
    tmax1 = tmax + et;
    xmin = xmin - ex;
    xmax = xmax + ex;
    ymin = ymin - ey;
    ymax = ymax + ey;
    
    axes(axy);
    switch type
        case 1
            plot(t,x,'color',color.tx,'linestyle','-');
            view(2);
            axis([tmin,tmax,xmin,xmax])
            xlabel('t')
            ylabel(xstring)
            set(r3dh,'vis','off')
            set(fig,'WindowButtonDownFcn','matpplane(''plxybdf'')');
        case 2
            plot(t,y,'color',color.ty,'linestyle','-');
            view(2);
            axis([tmin,tmax,ymin,ymax])
            xlabel('t')
            ylabel(ystring)
            set(r3dh,'vis','off')
            rotate3d off
            set(fig,'WindowButtonDownFcn','matpplane(''plxybdf'')');
        case 3
            plot(t,x,'color',color.tx,'linestyle','-');
            plot(t,y,'color',color.ty,'linestyle','--');
            axis([tmin,tmax,min(xmin,ymin),max(xmax,ymax)])
            view(2);
            xlabel('t')
            ylabel([xstring,' and ',ystring])
            set(r3dh,'vis','off')
            rotate3d off
            set(fig,'WindowButtonDownFcn','matpplane(''plxybdf'')');
        case 4
            plot3(x,y,t,'color',color.tx,'linestyle','-');
            view([-30,20]);
            axis([xmin,xmax,ymin,ymax,tmin1,tmax1]);
            zlabel('t')
            xlabel(xstring);
            ylabel(ystring);
            set(r3dh,'vis','on')
            rotate3d off
            set(fig,'WindowButtonDownFcn','');
        otherwise
            ax = [xmin,xmax,ymin,ymax,tmin1,tmax1];
            view([-30,20]);
            axis(ax);
            plot3(x,y,ax(5)*ones(size(x)),'color',color.ty,'linestyle','-');
            plot3(x,ax(4)*ones(size(y)),t,'color',color.ty,'linestyle','-');
            plot3(ax(2)*ones(size(x)),y,t,'color',color.ty,'linestyle','-');
            plot3(x,y,t,'color',color.tx,'linestyle','-');
            
            zlabel('t')
            xlabel(xstring);
            ylabel(ystring);
            set(r3dh,'vis','on')
            rotate3d off
            set(fig,'WindowButtonDownFcn','');
    end
    if (type < 4)
        set(aud.crop,'vis','on','enable','off');
    else
        set(aud.crop,'vis','off');
    end
    if type < 4
        set(fig,'windowbuttondownfcn','matpplane(''plxybdf'')');
    elseif type == 5
        set(fig,'windowbuttondownfcn',' ');
    end
    
elseif strcmp(action,'crop')
    
    cb = gcbo;
    axy = get(cb,'UserData');
    aud = get(axy,'UserData');
    delete(aud.h);
    set(aud.crop,'enable','off');
    aud.h = [];
    set(axy,'UserData',aud);
    ppxy = gcf;
    ud = get(ppxy,'UserData');
    data = ud.data;
    %color = ud.color;
    t = data(1,:);
    x = data(2,:);
    y = data(3,:);
    int = aud.int;
    k = find((t>=int(1)) & (t<=int(2)));
    t = t(k);
    x = x(k);
    y = y(k);
    type = get(aud.rad,'max');
    hh = matpplane('plotxyfig',type,ppxy);
    ud = get(hh(1),'UserData');
    ud.data = [t;x;y];
    set(hh(1),'UserData',ud);
    aud = get(hh(2),'UserData');
    matpplane('plxy',aud.rad,hh(1))
    
elseif strcmp(action,'plxybdf')
    
    dispa = gca;
    dispf = gcf;
    set(dispf,'windowbuttonmotionfcn','matpplane(''plxybmf'')',...
        'windowbuttonupfcn','matpplane(''plxybuf'')',...
        'Interrupt','on');
    aud = get(dispa,'UserData');
    delete(aud.h);
    aud.h = [];
    set(aud.crop,'enable','off');
    point = get(dispa,'currentpoint');
    aud.start = point(1,1);
    aud.finish = point(1,1);
    %  aud.centh = plot(point(1,1),point(1,2),'or','erase','xor');
    aud.h = plot(point(1,1),point(1,2),'--g',...
        'erase','xor','vis','off');
    set(dispa,'UserData',aud);
    
elseif strcmp(action,'plxybmf')
    
    dispa = gca;
    aud = get(dispa,'UserData');
    point = get(dispa,'currentpoint');
    finish = point(1,1);
    start = aud.start;
    xlim = get(dispa,'xlim');
    ylim = get(dispa,'ylim');
    
    if abs(finish - start)>0.05*(xlim(2)-xlim(1))
        set(aud.h,'xdata',[start,start,NaN,finish,finish],...
            'ydata',[ylim,NaN,ylim],...
            'vis','on');
    end
    aud.finish = finish;
    set(gca,'UserData',aud);
    
elseif strcmp(action,'plxybuf')
    
    dispa = gca;
    %dispf = gcf;
    xlim = get(dispa,'xlim');
    %ylim = get(dispa,'ylim');
    aud = get(dispa,'UserData');
    set(gcf,'windowbuttonmotionfcn','',...
        'windowbuttonupfcn','',...
        'Interrupt','on');
    start = aud.start;
    finish = aud.finish;
    if abs(finish - start)>0.05*(xlim(2)-xlim(1))
        set(aud.h,'erase','normal','vis','on');
        set(aud.crop,'enable','on');
        aud.int = [min(start,finish),max(start,finish)];
    else
        set(aud.h,'vis','off');
    end
    set(gca,'UserData',aud);
    
    
elseif strcmp(action,'print')
    
    dud = get(gcf,'UserData');
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'Preparing to print the matpplane Display Window.  Please be patient.';
    set(dud.notice,'string',nstr);
    
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'Printing the matpplane Display Window.';
    set(dud.notice,'string',nstr);
    set(gcf,'pointer','watch');
    eval(dud.printstr);
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'Ready.';
    set(dud.notice,'string',nstr);
    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.             ';];
    dud = get(gcf,'UserData');
    EqPtList = dud.eqpts;
    if isempty(EqPtList)
        disp('No equilibrium points have been computed.'),disp(' ')
    else
        disp('The following equilibrium points have been calculated:')
        disp(' ')
        L = size(EqPtList,1);
        %eqpttext = cell(L,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')
    
    disph = gcf;
    dud = get(disph,'UserData');
    axh = dud.axes;
    Xname = dud.syst.xvar;
    Yname = dud.syst.yvar;
    wmat = dud.wmat;
    WINvect = dud.syst.wind;
    
    NN = size(wmat,1);
    
    wch = 0;j=0;
    while wch == 0
        j = j+1;
        if WINvect == wmat(j,:)
            wch = j;
        end
    end
    winstr = cell(1,NN);
    for j = 1:NN
        a = num2str(wmat(j,1));
        b = num2str(wmat(j,2));
        c = num2str(wmat(j,3));
        d = num2str(wmat(j,4));
        winstr{j} = [' ',a,' < ',Xname,' < ',b,'   &   ',c,' < ', Yname,' < ',d];
    end
    
    [sel,~] = listdlg('liststring',winstr,...
        'selectionmode','single',...
        'listsize',[400,200],...
        'initialvalue',wch,...
        'name','matpplane Zoomback',...
        'promptstring','Select a rectangle:',...
        'OKString','Zoom');
    
    if (~isempty(sel))
        WINvect = wmat(sel,:);
        dud.syst.wind = WINvect;
        set(gcf,'UserData',dud);
        ppset = findobj('name','matpplane Setup');
        sud = get(ppset,'UserData');
        set(disph,'UserData',dud);
        set(sud.h.wind(1),'string',num2str(WINvect(1)));
        set(sud.h.wind(2),'string',num2str(WINvect(2)));
        set(sud.h.wind(3),'string',num2str(WINvect(3)));
        set(sud.h.wind(4),'string',num2str(WINvect(4)));
        sud.c.wind = WINvect;
        sud.o.wind = WINvect;
        set(ppset,'UserData',sud);
        aud = get(axh,'UserData');
        aud.DY = [WINvect(2) - WINvect(1);WINvect(4) - WINvect(3)];
        dwind = [WINvect(1); WINvect(3); -WINvect(2); -WINvect(4)];
        aud.cwind = dwind - dud.settings.magn*[aud.DY;aud.DY];
        set(axh,'UserData',aud);
        matpplane('dirfield',disph);
    end
    
elseif strcmp(action,'figdefault')
    % Opens the setup window. Initially called by "initialize" but called
    % many other places too
    fig = input1;
    set(fig,'CloseRequestFcn','matpplane(''closefcn'')');
    ppset = findobj('name','matpplane Setup');
    sud = get(ppset,'UserData');
    ud = get(fig,'UserData');
    ud.ssize = sud.ssize;
    fs = sud.fontsize;
    ud.fontsize = fs;
    style = sud.style;
    set(fig,'defaulttextfontsize',fs);
    set(fig,'defaultaxesfontsize',fs);
    set(fig,'defaultuicontrolfontsize',9*fs/10)
    lw = 0.5*fs/10;
    set(fig,'defaultaxeslinewidth',lw)
    set(fig,'defaultlinelinewidth',lw)
    set(fig,'defaultaxesfontname','helvetica')
    set(fig,'defaultaxesfontweight','normal')
    
    switch style % Note that an earlier line explicitly sets style = 'white' so none of this gets called
        case 'black'
            % if isunix | isvms, gamma = 0.5; else gamma = 0.0; end
            %whitebg(fig,[0,0,0])
            if isunix || isvms
                fc = [.35 .35 .35];
            else
                fc = [.2 .2 .2];
            end
            set(fig,'color',fc);
            set(fig,'defaultaxescolor',[0 0 0])
            % %whitebg(fig,brighten([.2 .2 .2],gamma))
            set(fig,'defaultaxescolor',[0 0 0])
            set(fig,'defaultaxescolororder', ...
                1-[0 0 1;0 1 0;1 0 0;0 1 1;1 0 1;1 1 0;.25 .25 .25]) % ymcrgbw
            % set(fig,'colormap',brighten(jet(64),gamma))
            set(fig,'colormap',jet(64))
            set(fig,'defaultsurfaceedgecolor',[0 0 0]);
        case 'white'
            %whitebg(fig,[1 1 1])
            set(fig,'color',[.8 .8 .8])
            set(fig,'defaultaxescolor',[1 1 1])
            set(fig,'defaultaxescolororder', ...
                [0 0 1;0 .5 0;1 0 0;0 .75 .75;.75 0 .75;.75 .75 0;.25 .25 .25]) % bgrymck
            set(fig,'colormap',jet(64))
            set(fig,'defaultsurfaceedgecolor',[0 0 0])
            
        case 'display'
            %whitebg(fig,[1 1 1])
            set(fig,'defaultaxescolor',[1 1 1])
            set(fig,'defaultaxescolororder', ...
                [0 0 1;0 .5 0;1 0 0;0 .75 .75;.75 0 .75;.75 .75 0;.25 .25 .25]) % bgrymck
            set(fig,'colormap',jet(64))
            set(fig,'defaultsurfaceedgecolor',[0 0 0])
            set(fig,'color',[1 1 1]*240/255);
            set(fig,'defaultuicontrolbackgroundcolor',[1 1 1]*220/255);
            set(fig,'defaultaxesfontweight','bold')
            set(fig,'defaulttextfontweight','bold')
            set(fig,'defaultaxeslinewidth',1)
            set(fig,'defaultlinelinewidth',1)
            
        case 'bw'
            %whitebg(fig,[0 0 0])
            set(fig,'color',[0 0 0])
            set(fig,'defaultaxescolor',[0 0 0])
            set(fig,'defaultaxescolororder', ...
                [1 1 1])
            set(fig,'colormap',[1 1 1;0 0 0])
            set(fig,'defaultsurfaceedgecolor',[0 0 0])
            
    end
    set(fig,'UserData',ud);
    
elseif strcmp(action,'zoominsq')
    
    % 'zoominsq' is the CallBackback for the Zoom in square menu item.
    
    set(gcf,'WindowButtonDownFcn','matpplane(''zoomsqd'')',...
        'WindowButtonUpFcn','1;','Interrupt','on');
    set(gca,'Interrupt','on');
    dud = get(gcf,'UserData');
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = ['Pick a center and ',...
        'drag the mouse, or just click on a center.'];
    set(dud.notice,'string',nstr);
    
elseif strcmp(action,'zoomsqd')
    
    ppdispa = gca;
    aud = get(ppdispa,'UserData');
    point = get(ppdispa,'currentpoint');
    aud.center = point(1,[1,2]);
    aud.centh = plot(point(1,1),point(1,2),'or','erase','xor');
    aud.box = plot(point(1,1),point(1,2),'--r','erase','xor');
    set(ppdispa,'UserData',aud);
    set(gcf,'windowbuttonmotionfcn','matpplane(''zoomsqm'')',...
        'windowbuttonupfcn','matpplane(''zoomsqu'')');
    
    
elseif strcmp(action,'zoomsqm')
    
    ppdispa = gca;
    aud = get(ppdispa,'UserData');
    point = get(ppdispa,'currentpoint');
    point = point(1,[1,2])';
    un = get(ppdispa,'units');
    set(ppdispa,'units','pix');
    w = get(ppdispa,'pos'); w = w([3,4]);w = w(:);
    set(ppdispa,'units',un);
    cent = aud.center(:);
    lam = max(abs(point - cent)./w);
    v = lam*w;
    data = cent(:,[1 1 1 1 1]) + [v,[-v(1);v(2)], -v, [v(1);-v(2)],v];
    set(aud.box,'xdata',data(1,:),'ydata',data(2,:));
    
elseif strcmp(action,'zoomsqu')
    
    disph = gcf;
    dud = get(disph,'UserData');
    axh = dud.axes;
    aud = get(axh,'UserData');
    hand = [aud.centh;aud.box];
    set(hand,'erase','normal');
    delete(hand);
    DY = aud.DY(:);
    un = get(axh,'units');
    set(axh,'units','pix');
    w = dud.syst.wind;
    point = get(axh,'currentpoint');
    point = point(1,[1,2]);point = point(:);
    cent = aud.center(:);
    ww = get(axh,'pos'); ww = ww([3,4]); ww = ww(:);
    lamb = max(abs(point-cent)./ww);
    v = lamb*ww;
    if ~all(v > 0.01*DY)
        points = [w([1 1 2 2]);w([3 4 3 4])];
        lambs=zeros(1,4);
        for j=1:4
            lambs(j) = max(abs(points(:,j)-cent)./ww);
        end
        lamb = min(lambs);
        v = lamb*ww/4;
    end
    p1 = cent + v; p2 = cent - v;
    a = [p2';p1'];
    a = [min(a);max(a)];
    DY = (a(2,:) - a(1,:))';
    WINvect = a(:)';
    dud.syst.wind = WINvect;
    aud.DY = DY;
    dwind = [WINvect(1); WINvect(3); -WINvect(2); -WINvect(4)];
    aud.cwind = dwind - dud.settings.magn*[aud.DY;aud.DY];
    set(axh,'units',un,'UserData',aud);
    set(disph,'UserData',dud);
    set(disph,'WindowButtonDownFcn','matpplane(''down'')',...
        'WindowButtonMotionFcn','matpplane(''cdisp'')',...
        'WindowButtonUpFcn','');
    matpplane('dirfield',disph);
    ppset = findobj('name','matpplane Setup');
    if isempty(ppset)
        matpplane('confused');
    else
        sud = get(ppset,'UserData');
        sud.c.wind = WINvect;
        sud.o.wind = WINvect;
        set(sud.h.wind(1),'string',num2str(WINvect(1)));
        set(sud.h.wind(2),'string',num2str(WINvect(2)));
        set(sud.h.wind(3),'string',num2str(WINvect(3)));
        set(sud.h.wind(4),'string',num2str(WINvect(4)));
        set(ppset,'UserData',sud);
    end
    
elseif strcmp(action,'level')
    
    ppdisp = gcf;
    dud = get(ppdisp,'UserData');
    lfcn = dud.level;
    %     Xname = dud.syst.xvar;
    %     Yname = dud.syst.yvar;
    %     Xname(find(abs(Xname)==92))=[];  % Remove \s if any.
    %     Yname(find(abs(Yname)==92))=[];
    Xname = strrep(dud.syst.xvar,'\','');
    Yname = strrep(dud.syst.yvar,'\','');
    pplevel = findobj('name','matpplane Level sets');
    if ~isempty(pplevel)
        delete(pplevel)
    end
    pplevel = figure('name','matpplane Level sets',...
        'vis','off',...
        'numb','off','tag','matpplane');
    
    matpplane('figdefault',pplevel);
    set(pplevel,'menubar','none');
    
    lev.fr1 = uicontrol('style','frame');
    lev.fr2 = uicontrol('style','frame');
    
    inst1str = ['Enter the function in terms of the variables ',...
        Xname, ' and ', Yname,':'];
    
    lev.inst1 = uicontrol('style','text','horiz','left',...
        'string',inst1str);
    
    lev.lfcn = uicontrol('style','edit','horiz','center',...
        'string',lfcn,'CallBack','',...
        'background','w');
    
    lev.ch(3) = uicontrol('style','radio','horiz','center',...
        'min',0,'max',3,...
        'value',0,...
        'vis','on',...
        'string','Let matpplane decide.');
    
    lev.ch(2) = uicontrol('style','radio','horiz','center',...
        'min',0,'max',2,...
        'value',0,...
        'string','Select a point in the Display Window.');
    
    lev.inst2 = uicontrol('style','text','horiz','left',...
        'string',['Choose one of the following ways to',...
        ' choose level value(s):']);
    
    lev.ch(1)  = uicontrol('style','radio','horiz','center',...
        'min',0,'max',1,...
        'value',0,...
        'string','Enter a vector of level values.');
    
    lev.rhs = uicontrol('style','edit','horiz','center',...
        'string',' ','CallBack','');
    
    lev.proc = uicontrol('style','push',...
        'string','Proceed',...
        'CallBack','matpplane(''levcomp'')');
    
    lev.close = uicontrol('style','push',...
        'string','Close',...
        'CallBack','close');
    
    for i=1:3
        set(lev.ch(i),'UserData',lev.ch(:,[1:(i-1),(i+1):3]));
    end
    
    CallBackrad = [
        'me = get(gcf,''currentobject'');',...
        'kk = get(me,''max'');',...
        'col = get(me,''backg'');',...
        'set(get(me,''UserData''),''value'',0),',...
        'set(me,''value'',kk);',...
        'ud = get(gcf,''UserData'');',...
        'if kk == 1,',...
        '   set(ud.rhs,''enable'',''on'',''backg'',''w'');',...
        'else,',...
        '   set(ud.rhs,''enable'',''off'',''backg'',col);',...
        'end,'];
    
    set(lev.ch,'CallBack',CallBackrad);
    
    left = 2; varl = 300; %buttw = 60;
    nudge = 3;
    tab = 15;
    lines1 = 5;
    lines2 = 2;
    xex = get(lev.inst1,'extent');
    ht = xex(4)+nudge;
    frw = varl + 2*tab+ 2*nudge;
    fr1bot = 2*left + ht;
    fr1ht = lines1*(nudge + ht) + nudge;
    fr2bot = fr1bot + fr1ht;
    fr2ht = lines2*(nudge + ht) + nudge;
    vbot = fr1bot + nudge;
    ch1bot = vbot + nudge + ht;
    ch2bot = ch1bot + nudge + ht;
    ch3bot = ch2bot + nudge + ht;
    inst2bot = ch3bot + nudge + ht;
    fbot =  fr2bot + nudge;
    inst1bot = fbot + nudge + ht;
    fleft = left + nudge + tab;
    instleft = left + nudge;
    chleft = instleft + tab;
    vleft = chleft + tab;
    vw = (frw - 2*vleft);
    fr1wind = [left,fr1bot,frw,fr1ht];
    fr2wind = [left,fr2bot,frw,fr2ht];
    inst1wind = [instleft,inst1bot,varl,ht];
    inst2wind = [instleft,inst2bot,varl,ht];
    fwind = [fleft,fbot,varl,ht];
    ch1wind = [chleft,ch1bot,varl,ht];
    ch2wind = [chleft,ch2bot,varl,ht];
    ch3wind = [chleft,ch3bot,varl,ht];
    vwind = [vleft,vbot,vw,ht];
    figw = 2*left + frw;
    fight = 3*left + ht + fr1ht + fr2ht;
    figwind = [40, 300, figw, fight];
    buttw = frw/2;
    sep = (figw - 2*buttw)/3;
    closel = sep;
    procl = 2*sep+buttw;
    clwind = [closel,left,buttw,ht];
    procwind = [procl,left,buttw,ht];
    set(pplevel,'pos',figwind);
    set(lev.fr1,'pos',fr1wind);
    set(lev.fr2,'pos',fr2wind);
    set(lev.inst1,'pos',inst1wind);
    set(lev.inst2,'pos',inst2wind);
    set(lev.ch(1),'pos',ch1wind);
    set(lev.ch(2),'pos',ch2wind);
    set(lev.ch(3),'pos',ch3wind);
    set(lev.rhs,'pos',vwind);
    set(lev.lfcn,'pos',fwind);
    set(lev.proc,'pos',procwind);
    set(lev.close,'pos',clwind);
    set(lev.ch(3),'value',3);
    set(lev.rhs,'enable','off');
    
    child = get(pplevel,'children');
    set(pplevel,'vis','on','UserData',lev);
    set(child,'vis','on');
    
elseif strcmp(action,'levcomp')
    % Plot Level curves
    pplevel = gcf;
    ud = get(pplevel,'UserData');
    ppdisp = findobj('name','matpplane Display');
    dud = get(ppdisp,'UserData');
    ppset = findobj('name','matpplane Setup');
    sud = get(ppset,'UserData');
    ch = ud.ch;
    val = zeros(1,3);
    for kk = 1:3
        val(kk) = get(ch(kk),'value');
    end
    KK = max(val);
    lfcn = get(ud.lfcn,'string');
    l=length(lfcn);
    for  k = fliplr(strfind(lfcn,'.'))
        if (find('*/^' == lfcn(k+1)))
            lfcn = [lfcn(1:k-1), lfcn(k+1:l)];
        end
        l=l-1;
    end
    pnameh = sud.h.pname;
    pvalh = sud.h.pval;
    %pflag = zeros(1,4);
    %perr = [];
    %lfcn(find(abs(lfcn)==32))=[];
    lfcn = strrep(lfcn,' ','');
    for kk = 1:4
        pn = get(pnameh(kk),'string');
        pv = get(pvalh(kk),'string');
        if ~isempty(pn)
            %pn(find(abs(pn)==92))=[];
            pn = strrep(pn,'\','');
            if isempty(pv)
                %perr = pvalh(kk);
            else
                pv = strrep(pv,' ','');
                %pv(find(abs(pv)==32))=[];
                lfcn = matpplane('paraeval',pn,pv,lfcn);
            end
        end
    end
    l = length(lfcn);
    for k=fliplr(find((lfcn=='^')|(lfcn=='*')|(lfcn=='/')))
        lfcn = [lfcn(1:k-1) '.' lfcn(k:l)];
        l = l+1;
    end
    WINvect = dud.syst.wind;
    XxXxXx = WINvect(1) + rand*(WINvect(2)-WINvect(1));
    YyYyYy = WINvect(3) + rand*(WINvect(4)-WINvect(3));
    %     Xname = dud.syst.xvar;
    %     Yname = dud.syst.yvar;
    %     Xname(find(abs(Xname)==92))=[];  % Remove \s if any.
    %     Yname(find(abs(Yname)==92))=[];
    Xname = strrep(dud.syst.xvar,'\','');
    Yname = strrep(dud.syst.yvar,'\','');
    
    if ~isvarname(Xname)
        errmsg = [Xname ' is not a valid variable name in MATLAB'];
        fprintf('\a')
        errordlg(errmsg,'MATPPLANE error','on');
        return;
    end
   
    if ~isvarname(Yname)
        errmsg = [Yname ' is not a valid variable name in MATLAB'];
        fprintf('\a')
        errordlg(errmsg,'MATPPLANE error','on');
        return;
    end
    
    
%     err = 0;res = 1;
%     eval([Xname,'=XxXxXx;'],'err = 1;');
%     eval([Yname,'=YyYyYy;'],'err = 1;');
%     eval(['res = ',lfcn, ';'],'err = 1;');
%     if err || isempty(res)
    lfcn=str2func(['@(' Xname ',' Yname ')' lfcn]);
    try
        lfcn(XxXxXx,YyYyYy);
    catch
        errmsg = 'The function does not evaluate correctly.';
        fprintf('\a')
        errordlg(errmsg,'MATPPLANE error','on');
        return;
    end

    Xmin = WINvect(1);
    Xmax = WINvect(2);
    Ymin = WINvect(3);
    Ymax = WINvect(4);
    N = 50; k = 4;
    deltax=(Xmax - Xmin)/(N-1);
    deltay=(Ymax - Ymin)/(N-1);
    XXXg=Xmin + deltax*(-k:N+k);
    YYYg=Ymin + deltay*(-k:N+k);
    
    [Xx,Yy]=meshgrid(XXXg,YYYg);
    Xxx=Xx(:);
    Yyy=Yy(:);
%    Ww = zeros(size(Xx(:)));
%    eval([Xname,'=Xxx;'],'err = 1;');
%    eval([Yname,'=Yyy;'],'err = 1;');
%    eval(['Ww = ',lfcn, ';'])
   Ww = lfcn(Xxx,Yyy); 

    KKK = 3; %# of significant figures.
    
    switch KK
        case 1  % vector input
            rhs = get(ud.rhs,'string');
            rhs = str2num(rhs); %#ok<ST2NM> %Takes vector data. can't be made str2double
            
        case 2  % mouse input
            figure(ppdisp);
            %     [XX,YY] = ppginput(1);
            [XX,YY] = ginput(1);
            figure(pplevel);
%            eval([Xname,'=XX;'],'err = 1;');
%            eval([Yname,'=YY;'],'err = 2;');
%            eval(['rhs = ',lfcn, ';'],'err = 3');
            rhs = lfcn(XX,YY);
            LL = ceil(log10(abs(rhs)));
            rhs = round(10^(KKK-LL)*rhs);
            rhs = 10^(LL-KKK)*rhs;
            
        case 3  % matpplane input
            %             Www = Ww;
            %             kkk = find(isnan(Www));
            %             Www(kkk) = [];
            %             kkk = find(imag(Www));
            %             Www(kkk) = [];
            Www = Ww(imag(Ww)==0);
            Www = rmmissing(Www);
            MM = max(Www);
            mm = min(Www);
            LL = ceil(log10(MM-mm));
            NN = 7;  % Number of curves
            rhs = mm+(1:NN).^2*(MM-mm)/NN^2;
            rhs = round(10^(KKK-LL)*rhs);
            rhs = 10^(LL-KKK)*rhs;
            
    end
    
    Ww = reshape(Ww,N+2*k+1,N+2*k+1);
    lrhs = length(rhs);
    if lrhs == 0
        return
    elseif lrhs == 1
        rhs = [rhs,rhs];
    end
    
    figure(ppdisp);
    [Cm,hcont] = contour(Xx,Yy,Ww,rhs,'--');
    hlabel = clabel(Cm,hcont);
    %  set(hlabel,'fontsize',dud.fontsize,...
    %	     'color',dud.color.level,...
    %	     'rotation',0);
    set(hlabel,'fontsize',dud.fontsize,...
        'color',[1,0,0],...
        'rotation',0);
    set(hcont,'visible','on',...
        'color',dud.color.level,...
        'linestyle','-.');
    dud.contours = [dud.contours ;hcont;hlabel];
    set(ppdisp,'UserData',dud);
    
elseif strcmp(action,'restart')
    
    ppset = findobj('name','matpplane Setup');
    oldfiles = dir('pptp*.m');
    for k = 1:length(oldfiles)
        fn = oldfiles(k).name;
        fid = fopen(fn,'r');
        fgetl(fid);  % Skip over 1st line
        fgetl(fid);  % Skip over 2nd line
        ll = fgetl(fid);
        fclose(fid);
        if strcmp(ll,'%% Created by matpplane')
            delete(fn)
        end
    end
    h = findobj('tag','matpplane');
    delete(setdiff(h,ppset));
    sud = get(ppset,'UserData');
    sud.flag = 0;
    set(ppset,'UserData',sud);
    figure(ppset)
    
elseif strcmp(action,'quit')
    % Called when the user chooses "quit".
    % Remove temporary directory from path
    ppset = findobj('name','matpplane Setup');
    sud = get(ppset,'UserData');
    if sud.remtd
        rmpath(tempdir);
    end
    % Delete temporary files
    oldfiles = dir([tempdir,'pptp*.m']);
    for k = 1:length(oldfiles)
        fn = [tempdir,oldfiles(k).name];
        fid = fopen(fn,'r');
        fgetl(fid);  % Skip over 1st line
        fgetl(fid);  % Skip over 2nd line
        ll = fgetl(fid);
        fclose(fid);
        if strcmp(ll,'%% Created by matpplane')
            delete(fn)
        end
    end
    % Close windows
    h = findobj('tag','matpplane');
    delete(h);
    
elseif strcmp(action,'closefcn')
    % This gets called when you click the close button in any window
    % Brings up the quit dialog in Setup or Display window
    fig = gcf;
    name = get(fig,'name');
    if strcmp(name,'matpplane Setup') || strcmp(name,'matpplane Display')
        quest = ['Closing this window will cause all matpplane ',...
            'windows to close, and matpplane will stop.  ',...
            'Do you want to quit matpplane?'];
        butt = questdlg(quest,'Quit matpplane?','Quit','Cancel','Quit');
        if strcmp(butt,'Quit')
            matpplane('quit');
        end
    elseif strcmp(name,'matpplane Linearization')
        dud = get(fig,'UserData');
        fcn = dud.function;
        if exist(fcn,'file')
            delete([fcn,'.m']);
        end
        delete(findobj('label',name));
        delete(fig);
    else
        delete(findobj('label',name));
        delete(fig);
    end
    
elseif strcmp(action,'confused')
    % Called in many places to indicate error and force restart
    tstring = 'matpplane is totally confused';
    qstring = {['You will have to restart matpplane from '...
        'the beginning in order to ',...
        'do anything new.  However, it might be possible '...
        'to save the current system ',...
        'or the gallery to make your restart easier, '...
        'or it may be possible to ',...
        'print out a figure, if the appropriate '...
        'figures are visible.  In such a case ',...
        'it would be best to do nothing now.'];
        'What do you want to do?'};
    bstr1 = 'Quit and restart matpplane.';
    bstr2 = 'Just quit matpplane.';
    bstr3 = 'Do nothing.';
    answer = questdlg(qstring,tstring,bstr1,bstr2,bstr3,bstr1);
    if strcmp(answer,bstr1)
        delete(findobj('tag','matpplane'));
        matpplane;return
    elseif strcmp(answer,bstr2)
        delete(findobj('tag','matpplane'));
        return
    else
        return
    end
    
elseif strcmp(action,'export')
    
    % export is the CallBackback for the Export solution data item in the
    % Options menu.
    
    disph = gcf;
    dud = get(disph,'UserData');
    arr = dud.arr;
    lv = get(arr.lines,'vis');
    av = get(arr.arrows,'vis');
    if ~isempty(arr.hx)
        nv = get(arr.hx(1),'vis');
    elseif ~isempty(arr.hy)
        nv = get(arr.hx(1),'vis');
    else
        nv = zeros(1,0);
    end
    if ~isempty(arr.barrows)
        bv = get(arr.barrows,'vis');
    else
        bv = zeros(1,0);
    end
    handles = [arr.lines;arr.arrows;arr.hx;arr.hy;arr.barrows];
    
    set(handles,'vis','off');
    %oldCallBack = get(disph,'WindowButtonDownFcn');
    set(disph,'WindowButtonDownFcn','');
    trjh = dud.solhand;
    notice = dud.notice;
    switch length(trjh)
        case 0
            if ishghandle(notice)
                nstr = get(notice,'string');
                nstr(1:3) = nstr(3:5);
                nstr{4} = 'There are no solutions.';
                nstr{5} = 'Ready.';
                set(notice,'string',nstr);
            end
            th = [];
            
        case 1
            th = trjh;
        otherwise
            if ishghandle(notice)
                nstr = get(notice,'string');
                nstr(1:4) = nstr(2:5);
                nstr{5} = 'Select a solution with the mouse.';
                set(notice,'string',nstr);
            end
            %    ppginput(1);
            ginput(1);
            th = get(disph,'currentobject');
    end
    if isempty(th)
        if ishghandle(notice)
            nstr = get(notice,'string');
            nstr(1:3) = nstr(3:5);
            nstr{4} = 'The item selected is not a solution.';
            nstr{5} = 'Ready.';
            set(notice,'string',nstr);
        end
    else
        vars = evalin('base','who');
        no = 1;
        kk = 0;
        while no
            kk = kk + 1;
            vstr = ['ppdata',num2str(kk)];
            if ~any(strcmp(vars,vstr))
                no = 0;
            end
        end
        yname = dud.syst.yvar;
        if abs(yname(1)) == 92
            yname = yname(2:length(yname));
        end
        xname = dud.syst.xvar;
        if abs(xname(1)) == 92
            xname = xname(2:length(xname));
        end
        tname = 't';
        tval = get(th,'zdata');
        tval = tval(:);
        xval = get(th,'xdata');
        xval = xval(:);
        yval = get(th,'ydata');
        yval = yval(:);
        ivstr = struct(tname,tval,xname,xval,yname,yval);
        assignin('base',vstr,ivstr);
        if ishghandle(notice)
            nstr = get(notice,'string');
            nstr(1:3) = nstr(3:5);
            nstr{4} = ['The data has been exported as the structure ',...
                vstr,' with fields ', tname, ', ', xname, ', and ', yname,'.'];
            nstr{5} = 'Ready.';
            set(notice,'string',nstr);
        end
        
        
    end
    
    
    
    set(arr.lines,'vis',lv);
    set(arr.arrows,'vis',av);
    set([arr.hx;arr.hy],'vis',nv);
    set(arr.barrows,'vis',bv);
    set(disph,'UserData',dud);
    set(disph,'WindowButtonDownFcn','matpplane(''down'')');
    
    
elseif strcmp(action,'cdisp')
    % Displays the cursor's current position at the bottom of figure
    [~,ppdisp] = gcbo;
    dud = get(ppdisp,'UserData');
    cp = get(ppdisp,'currentpoint');
    fpos = get(ppdisp,'pos');
    ppax = dud.axes;
    xd = get(ppax,'xlim');
    yd = get(ppax,'ylim');
    apos = get(ppax,'pos');
    xp = xd(1) + (cp(1) - apos(1)*fpos(3))*(xd(2)-xd(1))/(apos(3)*fpos(3));
    yp = yd(1) + (cp(2) - apos(2)*fpos(4))*(yd(2)-yd(1))/(apos(4)*fpos(4));
    str = ['(',num2str(xp,3),', ',num2str(yp,3),')'];
    set(dud.ccwind,'string',str);
    
elseif strcmp(action,'periodic')
    
    % Find a periodic orbit.
    
    ppdisp = findobj('name','matpplane Display');
    dud = get(ppdisp,'UserData');
    dfcn = dud.function;
    direction = dud.dir;
    notice = dud.notice;
    settings = dud.settings;
    refine = settings.refine;
    tol = settings.tol;
    AA = -1e6;
    BB = 1e6;
    switch direction
        case 0
            intplus = [0, BB];
            intminus = [0, AA];
        case -1
            intplus = [0, 0];
            intminus = [0, AA];
        case 1
            intplus = [0, BB];
            intminus = [0, 0];
    end
    ppdispa = dud.axes;
    aud = get(ppdispa,'UserData');
    DY = aud.DY;
    aud.plot = 0;
    atol = tol*DY*1e-4';
    set(ppdispa,'UserData',aud);
    nstr = get(dud.notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'Choose a starting point with the mouse.';
    set(dud.notice,'string',nstr);
    %  z00 = ppginput(1);
    z00 = ginput(1);
    z0 = z00;
    %ptstr = [' (',num2str(z0(1),2), ', ', num2str(z0(2),2), ')'];
    solver = dud.settings.solver;
    opt = odeset('OutputFcn',@matppout,'Refine',refine,...
        'RelTol',tol,'Abstol',atol);
    switch solver
        case 'Dormand Prince'
            solh = @matppdp45;
            opt = ppdisp;
        case 'Runge-Kutta 4'
            solh = @matpprk4;
            opt = ppdisp;
        case 'ode45'
            solh = @ode45;
        case 'ode23'
            solh = @ode23;
        case 'ode113'
            solh = @ode113;
        case 'ode15s'
            solh = @ode15s;
        case 'ode23s'
            solh = @ode23s;
        case 'ode23t'
            solh = @ode23t;
        case 'ode23tb'
            solh = @ode23tb;
    end
    %exist(dfcn);  % Does nothing
    dfh = str2func(dfcn);
    dud.noticeflag = 0;   % Notices only from here.
    set(ppdisp,'UserData',dud);
    if intplus(2)>intplus(1)
        if ishghandle(notice)
            nstr = get(notice,'string');
            nstr(1:4) = nstr(2:5);
            nstr{5} = 'The forward orbit ';
            set(notice,'string',nstr);
        end
        drawnow
        [tp,xp] = solh(dfh,intplus,z0,opt);
        aud = get(ppdispa,'UserData');
        set(aud.line,'erase','normal')
        delete(aud.line);
        aud.line = [];
        z0 = aud.y;
        stopcode = aud.stop;
        nstr = get(dud.notice,'string');
        switch stopcode
            case 1
                nstr{5} = [nstr{5}, ' left the computation window.'];
            case 2
                zz = aud.zz;
                ystr = ['(',num2str(zz(1),2), ', ', num2str(zz(2),2),').'];
                nstr{5} = [nstr{5}, ' --> a possible eq. pt. near ',ystr];
            case 3
                nstr{5} = [nstr{5}, ' --> a closed orbit'];
            case 4
                nstr{5} = [nstr{5}, ' was stopped by the UserData.'];
            case 5
                y = aud.y; % RG: y was undefined. This seems like the right fix
                ystr = ['(',num2str(y(1),2), ', ', num2str(y(2),2),').'];
                nstr(1:3) = nstr(2:4);
                nstr{4} = [nstr{5},' experienced a failure at ',ystr];
                nstr{5} = 'Problem is singular or tolerances are too large.';
        end
        set(dud.notice,'string',nstr);
        drawnow
        
        if stopcode==3  % periodic orbit found.
            % [tp,xp] = feval(solh,dfh,intplus,z0,opt);
            rr = sqrt(((xp(:,1) - z0(1))/DY(1)).^2 + ((xp(:,2) - z0(2))/DY(2)).^2);
            kk = find(rr < 0.01);
            %LL = length(kk);
            kmax = max(kk);
            kkk = max(kk(kk<kmax-30));
            T = tp(kmax) - tp(kkk);
            aud.gstop = 0;
            set(ppdispa,'UserData',aud);
            kkk = 0;
            NN = 1;
            dz0 = feval(dfh,0,z0);
            while abs(NN)>0.0001*abs(T) && kkk < 10
                [tt,yy] = solh(dfh,[0,T],z0,opt);
                yy = yy(length(tt),:);
                dyy = feval(dfh,0,yy');
                NN = ((yy-z0')*dz0)/(dyy'*dz0);
                T = T - NN;
                kkk = kkk +1;
            end
            nstr{5} = [nstr{5}, ' of period ', num2str(abs(T),3), '.'];
            set(dud.notice,'string',nstr);
            drawnow
            [tp,xp] = solh(dfh,[0 3*T],z0,opt);
            aud.gstop = 1;
            set(ppdispa,'UserData',aud);
            hnew = plot(xp(:,1),xp(:,2),'color','b');
            set(hnew,'zdata',tp);
            solhand = [dud.solhand;hnew];
            dud.solhand = solhand;
        end
        drawnow
    end
    if intminus(2)<intminus(1)
        z0 = z00;
        if ishghandle(notice)
            nstr = get(notice,'string');
            nstr(1:4) = nstr(2:5);
            nstr{5} = 'The backward orbit ';
            set(notice,'string',nstr);
        end
        drawnow
        [tp,xp] = solh(dfh,intminus,z0,opt);
        aud = get(ppdispa,'UserData');
        set(aud.line,'erase','normal')
        delete(aud.line);
        aud.line = [];
        %z0 = aud.y;
        stopcode = aud.stop;
        nstr = get(dud.notice,'string');
        switch stopcode
            case 1
                nstr{5} = [nstr{5}, ' left the computation window.'];
            case 2
                zz = aud.zz;
                ystr = ['(',num2str(zz(1),2), ', ', num2str(zz(2),2),').'];
                nstr{5} = [nstr{5}, ' --> a possible eq. pt. near ',ystr];
            case 3
                nstr{5} = [nstr{5}, ' --> a nearly closed orbit.'];
            case 4
                nstr{5} = [nstr{5}, ' was stopped by the UserData.'];
            case 5
                ystr = ['(',num2str(y(1),2), ', ', num2str(y(2),2),').'];
                nstr(1:3) = nstr(2:4);
                nstr{4} = [nstr{5},' experienced a failure at ',ystr];
                nstr{5} = 'Problem is singular or tolerances are too large.';
        end
        set(dud.notice,'string',nstr);
        drawnow
        if stopcode==3  % periodic orbit found.
            z0 = aud.y;
            rr = sqrt(((xp(:,1) - z0(1))/DY(1)).^2 + ((xp(:,2) - z0(2))/DY(2)).^2);
            kk = find(rr < 0.01);
            %LL = length(kk);
            kmax = max(kk);
            kkk = max(kk(kk<kmax-30));
            T = tp(kmax) - tp(kkk);
            aud.gstop = 0;
            set(ppdispa,'UserData',aud);
            kkk = 0;
            NN = 1;
            dz0 = feval(dfh,0,z0);
            while abs(NN)>0.0001*abs(T) && kkk < 10
                [tt,yy] = solh(dfh,[0,T],z0,opt);
                yy = yy(length(tt),:);
                dyy = feval(dfh,0,yy');
                NN = ((yy-z0')*dz0)/(dyy'*dz0);
                T = T - NN;
                kkk = kkk +1;
            end
            nstr{5} = [nstr{5}, ' of period ', num2str(abs(T),3), '.'];
            set(dud.notice,'string',nstr);
            drawnow
            [tp,xp] = solh(dfh,[0 3*T],z0,opt);
            aud.gstop = 1;
            set(ppdispa,'UserData',aud);
            hnew = plot(xp(:,1),xp(:,2),'color','b');
            set(hnew,'zdata',tp);
            solhand = [dud.solhand;hnew];
            dud.solhand = solhand;
        end
        drawnow
    end
    dud.noticeflag = 1;
    dud.dir = 0;
    set(ppdisp,'UserData',dud);
    aud.plot = 1;
    set(ppdispa,'UserData',aud);
    
    if ishghandle(notice)
        nstr = get(notice,'string');
        nstr(1:4) = nstr(2:5);
        nstr{5} = 'Ready';
        set(notice,'string',nstr);
    end
    drawnow
end