function output = odesolve(action,input1)

% ODESOLVE is a MATLAB function which provides a Graphical User
%          Interface (GUI) for the use of MATLAB's differential
%          equation solvers.   

% Copyright (c) 2003 John C. Polking, Rice University.  

startstr = 'odesolve';

if nargin < 1
    action ='initialize';
end


switch action
    case 'initialize'
        
        % First we make sure that there is no other copy of ODESOLVE
        % running, since this causes problems.	
        
        dfh = findobj('tag','odesolve');	
        if ~isempty(dfh);
            qstring = {'There are some open ODESOLVE figures.  These may be' ...
                    ' invisible.   What do you want to do?'};
            tstring = 'Only one copy of ODESOLVE can be open at one time.';
            b1str = 'Close and restart.';
            b2str = 'Just close.';
            b3str = 'Nothing.';
            answer = questdlg(qstring,tstring,b1str,b2str,b3str,b1str);
            switch answer
                case b1str
                    delete(dfh);
                    eval(startstr);return
                case b2str
                    delete(dfh);return
                case b3str
                    return
            end 
        end  % if ~isempty(dfh);
        
        % 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;
        if isempty(findstr(tmpdir,p))
            ud.remtd = 1;
            addpath(tempdir)
        end
        
        % Next we look for old files created by odesolve.
        
        oldfiles = dir('odes*.m');
        kk = zeros(0,1);
        for k = 1:length(oldfiles)
            fn = oldfiles(k).name;
            fid = fopen(fn,'r');
            ll = fgetl(fid);
            ll = fgetl(fid);
            ll = fgetl(fid);
            fclose(fid);
            if strcmp(ll,'%% Created by odesolve')
                delete(fn)
            else
                kk = [kk;k];
            end
        end
        oldfiles = dir([tempdir,'odes*.m']);
        for k = 1:length(oldfiles)
            fn = [tempdir,oldfiles(k).name];
            fid = fopen(fn,'r');
            ll = fgetl(fid);
            ll = fgetl(fid);
            ll = fgetl(fid);
            fclose(fid);
            if strcmp(ll,'%% Created by odesolve')
                delete(fn)
            else
                kk = [kk;k];
            end
        end
        
        %%  Default settings.
        
        ud.ssize = get(0,'defaultaxesfontsize');
        ud.oddir = pwd;
        ud.solvers = {'ode45';
            'ode23';
            'ode113';
            'ode15s';
            'ode23s';
            'ode23t';
            'ode23tb'};
        ud.outputhandle = @tpout;
        comp = computer;
        if strcmp(comp,'PCWIN')
            ud.fontsize = 0.8*ud.ssize;
        else
            ud.fontsize = ud.ssize;
        end
        
        ecolor = 'w';  %% Background for text boxes
        dcolor = 'r';  %% Indicates errors
        ccolor = 'y';  %% Caution
        
        %% The initial window
        
        odesolveset = figure('name','ODESOLVE Setup','numb','off',...
            'tag','odesolve','visible','off',...
            'user',ud);
        tcolor = get(gcf,'defaultuicontrolbackgroundcolor');
        odesolve('figdefault',odesolveset)
        
        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;      % 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 = 1;		% Left margin of the frames.
        frsep = 1;    	% Separation between frames.
        separation = texth;	% Separation between boxes.
        eqleft = left +5;
        fudge = 0.15*separation;
        
        system.name = 'default system';
        system.NN = 2;
        system.xname = {'x', 'y'};
        system.der = {' y', ' y^2 - x' };
        system.ic = {'2', '0'};
        system.tname = 't';
        system.tvals = [0 0 30];
        system.pname = {};
        system.pval = [];
        system.exname = {};
        system.exval = {};
        
        system(2).name = 'forced oscillator';
        system(2).NN = 2;
        system(2).xname = {'x', 'v'};
        system(2).der = {' v', ' -C*v -om0^2*x + A*cos(om *t)' };
        system(2).ic = {'0', '0'};
        system(2).tname = 't';
        system(2).tvals = [0 0 30];
        system(2).pname = {'C', 'om0','A','om'};
        system(2).pval = [1,4,1,2];
        system(2).exname = {};
        system(2).exval = {};
        
        system(3).name = 'immune system';
        system(3).NN = 3;
        system(3).xname = {'E_1'  'E_2'  'V'};
        system(3).der = {' GG - mu*E_1 + EEEE + KK*V*E_1'  ' GG - mu*E_2 + EEEE'  ' rr*V - kk*V*E_1'};
        system(3).ic = {'1', '1', '1'};
        system(3).tname = 't';
        system(3).tvals = [0 0 100];
        system(3).pname = {'GG'  'mu'  'rr'  'KK'  'kk'  ''};
        system(3).pval = [1 1.25 0.1 0.05 0.01];
        system(3).exname = {'EEEE'  ''  ''  ''};
        system(3).exval = {'0.252*E_1*E_2/(1+0.008*E_1*E_2)'};
        
        system(4).name = 'forced Duffing';
        system(4).NN = 2;
        system(4).xname = {'x','v'};
        system(4).der = {' v',' (-C*v -om0^2*x - k*x^3 + FF)/m'};
        system(4).ic = {'0', '0'};
        system(4).tname = 't';
        system(4).tvals = [0   0  30];
        system(4).pname = {'C','om0','k','om','A','m'};
        system(4).pval = [1  4  1  2  1 1];
        system(4).exname = {'FF','','',''};
        system(4).exval = {'A*cos(om *t)'};
        
        
        system(5).name = 'Lorenz';
        system(5).NN = 3;
        system(5).xname = {'x','y','z'};
        system(5).der = {'  -a*x + a*y',' b*x - y - x*z',' -r*z + x*y'};
        system(5).ic = {'10', '-10',  '15'};
        system(5).tname = 't';
        system(5).tvals = [0   0 100];
        system(5).pname = {'a','r','b','','',''};
        system(5).pval = [10      2.66667           28];
        system(5).exname = {'','','',''};
        system(5).exval = {};
        
        system(6).name = 'logistic';
        system(6).NN = 1;
        system(6).xname = {'P'};
        system(6).der = {'rr*P*(1 - P)'};
        system(6).ic = {'0.1'};
        system(6).tname = 't';
        system(6).tvals = [0   0 20];
        system(6).pname = {'rr'};
        system(6).pval = [0.5];
        system(6).exname = {'','','',''};
        system(6).exval = {};
        
        system(7).name = '1D stiff';
        system(7).NN = 1;
        system(7).xname = {'y'};
        system(7).der = {'exp(t)*cos(y)'};
        system(7).ic = {'0'};
        system(7).tname = 't';
        system(7).tvals = [0   0 12];
        system(7).pname = {};
        system(7).pval = [];
        system(7).exname = {'','','',''};
        system(7).exval = {};
        
        H.name = 'two springs';
        H.NN = 4;
        H.xname = {'x','u','y','v'};
        H.der = {'u','(k2*(y-x) - k1*x)/M1','v','(-k2*(y-x) + F)/M2'};
        H.ic = {'2','0','2','0'};
        H.tname = 't';
        H.tvals = [0   0  30];
        H.pname = {'k1','k2','A','M1','M2','om'};
        H.pval = [1  1  0  1  2  4];
        H.exname = {'F','','',''};
        H.exval = {'A*cos(om*t)','','',''};
        
        system(8) = H;
        
        ud.c = system(1);	% Changed values.
        ud.o = system(1);	% Original values.
        %%  ud.h.NN = system(1).NN;
        ud.h.name = system(1).name;
        ud.h.xname = [];  % This will be the handles in the
        ud.h.equals = [];
        ud.h.der = [];
        ud.h.icname = [];
        ud.h.icval = [];
        ud.h.pname = [];
        ud.system = system;
        
        eqfrw = varw+equalw+eqlength+10;
        
        icnamew = 38*rr;
        icvalw = 60*rr;
        % iceqw = 8*rr;
        % icfrw = icnamew+iceqw+icvalw+10;
        icfrw = icnamew+icvalw+10;
        setfigw = 2*left + eqfrw + icfrw + frsep;
        
        %%  The buttons
        
        buttw = setfigw/3;
        qwind = [0,frsep,buttw,separation];	        % Quit button
        rwind = [buttw,frsep,buttw,separation];	% Revert
        sowind = [buttw,frsep,buttw,separation];	% Solve in old fig
        snwind = [2*buttw,frsep,buttw,separation];	% Solve in new fig
        buttht = 2*frsep+separation;
        
        butt(1) = uicontrol('style','push',...
            'pos',qwind,...
            'string','Quit',...
            'call','odesolve(''quit'')',...
            'visible','off');
        
        socall = [
            'ud = get(gcf,''user'');'...
                'ud.addflag = 1;'...
                'set(gcf,''user'',ud);'...
                'odesolve(''solve'')'];
        
        butt(2) = uicontrol('style','push',...
            'pos',sowind,...
            'string','Add to the existing display',...
            'call',socall,...
            'enable','off',...
            'visible','off');
        ud.h.sobutt = butt(2);
        
        sncall = [
            'ud = get(gcf,''user'');'...
                'ud.addflag = 0;'...
                'set(gcf,''user'',ud);'...
                'odesolve(''solve'')'];
        
        butt(3) = uicontrol('style','push',...
            'pos',snwind,...
            'string','Solve in a new window',...
            'call',sncall,....
            'visible','off');
        
        %% The Options frame
        
        opfrbot = buttht;
        opfrw = setfigw - 2*left;
        opfrht = separation + 10;
        opfrwind = [left, opfrbot, opfrw, opfrht];
        opframeht = opfrbot + opfrht + frsep;
        opframe = uicontrol('style','frame','pos',opfrwind);
        opw = (setfigw - 2*eqleft)/5;
        opbot = opfrbot + 5;
        solvtith = uicontrol('style','text',...
            'pos',[eqleft,opbot-fudge,opw,texth],...
            'horizon','right',...
            'string','Solver:    ',...
            'visible','off',...
            'backgroundcolor',tcolor);
        solvcall = [
            '[h,fig] = gcbo;',...
                'sud = get(fig,''user'');',...
                'oldsolver = sud.solver;',...
                'solverlist = get(h,''string'');',...
                'newsolver = solverlist{get(h,''value'')};',...
                'if ~strcmp(oldsolver,newsolver),',...
                '  sud.settings = sud.defaultsettings;',...
                '  if ~strcmp(newsolver,''ode45''),',...
                '    sud.settings.refine = 1;',...
                '  end,',...
                '  sud.solver = newsolver;',...
                '  set(fig,''user'',sud);',...
                'end'];
        
        ud.h.solver = uicontrol('style','popup',...
            'pos',[eqleft+opw opbot, opw texth],...
            'horizon','center',...
            'string',ud.solvers,...
            'user','',...
            'call',solvcall,...
            'visible','off',...
            'backgroundcolor',ecolor);
        outtith = uicontrol('style','text',...
            'pos',[eqleft+2*opw,opbot-fudge,opw,texth],...
            'horizon','right',...
            'string','Output:    ',...
            'visible','off',...
            'backgroundcolor',tcolor);
        outputlist = {'time plot', '2D phase plot', '3D plot'};
        ud.h.output = uicontrol('style','popup',...
            'pos',[eqleft+3*opw opbot, opw texth],...
            'horizon','center',...
            'string',outputlist,...
            'user',1,...
            'call','odesolve(''outputchoice'')',...
            'visible','off',...
            'backgroundcolor',ecolor);
        ud.h.outputstr = uicontrol('style','text',...
            'pos',[left+4*opw,opbot-fudge,opw,texth],...
            'horizon','left',...
            'string','     all variables vs t',...
            'visible','off',...
            'backgroundcolor',tcolor);
        
        
        
        %% The expressions frame
        
        exfrbot = opframeht;
        exfrw = setfigw - 2*left;
        exfrht = 2*separation + 10;
        exfrwind = [left, exfrbot, exfrw, exfrht];
        exframeht = exfrbot + exfrht + frsep;
        exframe = uicontrol('style','frame','pos',exfrwind);
        exnamew = 30*rr;
        exeqw = 8*rr;
        exbot(2) = exfrbot + 5;
        exbot(1) = exbot(2) +separation;
        extit = uicontrol('style','text',...
            'horizon','center',...
            'string','Expressions:',...
            'visible','off','backgroundcolor',tcolor);
        ext = get(extit,'extent');
        extitw = ext(3);
        pos = [eqleft exbot(1)-fudge extitw texth];
        set(extit,'pos',pos);
        exncall = [
            '[h,fig] = gcbo;'...
                'ud =  get(fig,''user'');'...
                'num = get(h,''user'');'...
                'ud.c.exname{num} = get(ud.h.exname(num),''string'');'...
                'ud.fileflag = 1;'...
                'odesolve(''reset'');'...
                'set(fig,''user'',ud);'];
        exvcall = [
            '[h,fig] = gcbo;'...
                'ud =  get(fig,''user'');'...
                'num = get(h,''user'');'...
                'ud.c.exval{num} = get(ud.h.exval(num),''string'');'...
                'ud.fileflag = 1;'...
                'set(h,''backgroundcolor'',ud.ecolor);'...
                'odesolve(''reset'');'...
                'set(fig,''user'',ud);'];
        exsep = 20;
        exvalw = (setfigw - 2*eqleft - extitw)/2 - exsep - exnamew - exeqw;
        exval = ud.c.exval;
        exname = ud.c.exname;
        for jj = 1:2
            for kk = 1:2
                exleft = eqleft +extitw +exsep +(kk-1)*(exnamew+exeqw+exvalw+ exsep);
                exeqleft = exleft+exnamew;
                exvleft = exeqleft +exeqw;
                K = kk + 2*(jj-1);
                if K > length(exname)
                    exname{K} = '';
                    exval{K} = '';
                end
                name = exname{K};
                value = exval{K};
                ud.h.exname(K) = uicontrol('style','edit',...
                    'pos',[exleft exbot(jj) exnamew texth],...
                    'horizon','right','string',name,...
                    'user',K,...
                    'call',exncall,...
                    'visible','off',...
                    'backgroundcolor',ecolor);
                
                exequal(K) = uicontrol('style','text',...
                    'pos',[exeqleft exbot(jj)-fudge exeqw texth],...
                    'horizon','center',...
                    'string','=',...
                    'visible','off',...
                    'backgroundcolor',tcolor);
                
                ud.h.exval(K) = uicontrol('style','edit',...
                    'pos',[exvleft exbot(jj) exvalw texth],...
                    'string',value,...
                    'call',exvcall,...
                    'visible','off',...
                    'user',K,...
                    'backgroundcolor',ecolor);
            end
        end
        ud.c.exname = exname;
        ud.c.exval = exval;
        
        
        
        %% The parameter frame
        
        pfrbot = exframeht;
        pfrw = setfigw - 2*left;
        pfrht = 2*separation + 10;
        pfrwind = [left, pfrbot, pfrw, pfrht];
        pframeht = pfrbot + pfrht + frsep;
        pframe = uicontrol('style','frame','pos',pfrwind);
        pncall = [
            '[h,fig] = gcbo;'...
                'ud =  get(fig,''user'');'...
                'num = get(h,''user'');'...
                'ud.c.pname{num} = get(ud.h.pname(num),''string'');'...
                'ud.fileflag = 1;'...
                'odesolve(''reset'');'...
                'set(gcf,''user'',ud);'];
        
        pvcall = [
            '[h,fig] = gcbo;'...
                'ud =  get(fig,''user'');'...
                'num = get(h,''user'');'...
                'nstr = str2num(get(ud.h.pval(num),''string''));'...
                'if ~isempty(nstr),'...
                'ud.c.pval(num) = str2num(get(ud.h.pval(num),''string''));'...
                'else, ud.c.pval(num) = NaN;'...
                'end,'...
                'set(h,''backgroundcolor'',ud.ecolor);'...
                'set(gcf,''user'',ud);'];
        
        pnamew = 30*rr;
        peqw = 8*rr;
        pbot1 = pfrbot + 5;
        pbot2 = pbot1 +separation;
        pbot(1) = pbot2;
        pbot(2) = pbot1;
        
        paratit=uicontrol('style','text',...
            'horizon','center',...
            'string','Parameters:',...
            'visible','off','backgroundcolor',tcolor);
        
        paratitw = extitw;
        pos = [eqleft pbot2-fudge paratitw texth];
        set(paratit,'pos',pos);
        psep = 20;
        pvalw = (setfigw - 2*eqleft - paratitw)/3 - psep - pnamew - peqw;
        pval = ud.c.pval;
        pname = ud.c.pname;
        for jj = 1:2
            for kk = 1:3
                pleft = eqleft + paratitw +psep +(kk-1)*(pnamew+peqw+pvalw+ psep);
                peqleft = pleft + pnamew;
                pvleft = peqleft + peqw;
                K = kk +3*(jj-1);
                if K > length(pname)
                    pname{K} = '';
                end
                name = pname{K};
                if K <= length(pval)
                    value = num2str(pval(K));
                else
                    value = '';
                end
                ud.h.pname(K) = uicontrol('style','edit',...
                    'pos',[pleft pbot(jj) pnamew texth],...
                    'horizon','right','string',name,...
                    'user',K,...
                    'call',pncall,...
                    'visible','off',...
                    'backgroundcolor',ecolor);
                equal(K) = uicontrol('style','text',...
                    'pos',[peqleft pbot(jj)-fudge peqw texth],...
                    'horizon','center',...
                    'string','=',...
                    'visible','off',...
                    'backgroundcolor',tcolor);
                
                ud.h.pval(K) = uicontrol('style','edit',...
                    'pos',[pvleft pbot(jj) pvalw texth],...
                    'string',value,...
                    'call',pvcall,...
                    'visible','off',...
                    'user',K,...
                    'backgroundcolor',ecolor);
            end
        end
        ud.c.pname = pname;
        ud.c.pval = pval;
        
        %% The number of equations frame
        
        frw = (eqfrw-2)/3;  %% The width of the middle three frames.
        
        nfrbot = pframeht;
        nfrw = frw;
        nfrht = 2*separation + 10;
        nfrwind = [left, nfrbot, nfrw, nfrht];
        nframe = uicontrol('style','frame','pos',nfrwind);
        nbot1 = nfrbot + 5;
        nbot2 = nbot1 + separation;
        numbeqw = nfrw - 4;
        numbeqwind = [left+2 nbot2 numbeqw texth];
        numbeq = uicontrol('style','text',...
            'pos',numbeqwind,...
            'string','Number of equations.',...
            'background',tcolor);
        
        neqnw = 30*rr;
        ll = (numbeqw - neqnw)/2;
        neqnwind = [2+ll nbot1 neqnw texth];
        neqn = uicontrol('style','edit',...
            'pos',neqnwind,...
            'string','2',...
            'call','odesolve(''numchange'')',...
            'background',ecolor);
        ud.h.NN = neqn;
        
        %% The independent variable frame.
        
        tfrbot = pframeht;
        tfrleft = left + nfrw + frsep;
        tfrw = frw;
        tfrht = 2*separation + 10;
        tfrwind = [tfrleft,tfrbot,tfrw,tfrht];
        tframeht = tfrbot + tfrht + frsep;
        tframe = uicontrol('style','frame','pos',tfrwind);
        
        tbot1 = nbot1;
        tbot2 = nbot2;
        tstrw = tfrw - 4;
        tstrwind = [tfrleft+2, tbot2, tstrw, texth];
        tstr = uicontrol('style','text',...
            'pos',tstrwind,...
            'string','Independent variable.',...
            'background',tcolor);
        tvarw = 30*rr;
        ll = (tstrw - tvarw)/2;
        tvarwind = [tfrleft+2+ll tbot1 tvarw texth];
        tcall = [
            '[h,fig] = gcbo;'...
                'ud =  get(fig,''user'');'...
                'str = get(h,''string'');'...
                'ud.c.tname = str;'...
                'set(ud.h.tic,''string'',[str '' = '']);'...
                'set(ud.h.solint(1),''string'',['' <= '',str, '' <= '']);'...
                'str = [''     all variables vs '', str, ''.''];'...
                'set(ud.h.outputstr,''string'',str);'...
                'ud.fileflag = 1;'...
                'set(fig,''user'',ud);' ];
        
        
        ud.h.tvar = uicontrol('style','edit',...
            'pos',tvarwind,...
            'string','t',...
            'call',tcall,...
            'background',ecolor);
        
        
        
        %% The solution interval frame.
        
        solfrbot = pframeht;
        solfrleft = tfrleft + tfrw + frsep;
        solfrw = frw;
        solfrht = 2*separation + 10;
        solfrwind = [solfrleft,solfrbot,solfrw,solfrht];
        solframeht = tfrbot + tfrht + frsep;
        solframe = uicontrol('style','frame','pos',solfrwind);
        
        solbot1 = nbot1;
        solbot2 = nbot2;
        solstrw = solfrw - 4;
        solstrwind = [solfrleft+2, solbot2, solstrw, texth];
        solstr = uicontrol('style','text',...
            'pos',solstrwind,...
            'string','Solution interval.',...
            'background',tcolor);
        solvarw = 30*rr;
        ll = (solstrw - solvarw)/2;
        solvarwind = [solfrleft+2+ll tbot1 tvarw texth];
        
        solineqstr = ['  <=  ',ud.c.tname, '  <=  '];
        ud.h.solint(1) = uicontrol('style','text',...
            'string',solineqstr,...
            'background',tcolor);
        ext = get(ud.h.solint(1),'extent');
        www = ext(3);
        ddd = (solfrw - www)/2 ;
        solineql = ddd + solfrleft;
        solineqwind = [solineql, tbot1-fudge ,www,texth];
        set(ud.h.solint(1),'pos',solineqwind);
        wind = [ solfrleft+3 tbot1,ddd-3,texth];
        str = num2str(ud.c.tvals(2));
        solintcall = [
            '[h,fig] = gcbo;'...
                'ud =  get(fig,''user'');'...
                'tvals = ud.c.tvals;'...
                'num = get(h,''user'');'...
                'value = str2num(get(h,''string''));'...
                'tvals(num) = value;'...
                'if (tvals(2)>tvals(1) | tvals(3)<tvals(1)),',...
                '   col = ud.dcolor;',...
                'else,',...
                '   col = ud.ecolor; end,'...
                'set([ud.h.solint(2:3) ud.h.tval],''background'',col);'...
                'ud.c.tvals =tvals;'...
                'set(fig,''user'',ud);'];
        
        ud.h.solint(2) = uicontrol('style','edit',...
            'pos',wind,...
            'horizon','right',...
            'string',str,...
            'call',solintcall,...
            'user',2,...
            'visible','off',...
            'backgroundcolor',ecolor);
        wind = [ solineql+www tbot1,ddd-3,texth];
        str = num2str(ud.c.tvals(3));
        ud.h.solint(3) = uicontrol('style','edit',...
            'pos',wind,...
            'horizon','left',...
            'string',str,...
            'call',solintcall,...
            'user',3,...
            'visible','off',...
            'backgroundcolor',ecolor);
        
        %% The equation frame.
        
        eqfrbot = tframeht;
        eqfrht = 3*separation + 10;
        eqfrwind = [left, eqfrbot, eqfrw, eqfrht];
        frame(1) = uicontrol('style','frame','pos',eqfrwind);
        ud.eqnlabel = uicontrol('style','text',...
            'string','The differential equations.',...
            'vis','off',...
            'background',tcolor);
        
        %% The ic frame
        
        icfrbot = solfrbot;
        icfrht = solfrht + eqfrht +frsep;
        icfrl = 2*left+eqfrw;
        icnleft = icfrl +2;
        icvleft = icnleft + icnamew;
        icfrwind = [icfrl, icfrbot, icfrw, icfrht]; 
        frame(2) = uicontrol('style','frame','pos',icfrwind);
        ud.iclabel = uicontrol('style','text',...
            'string','The initial conditions.',...
            'vis','off',...
            'background',tcolor);
        tbot = nbot1;
        tnamstr = [ud.c.tname, ' = '];
        tvalstr = num2str(ud.c.tvals(1));
        ud.h.tic = uicontrol('style','text',...
            'pos',[icnleft,tbot-fudge,icnamew,texth],...
            'horizon','right',...
            'string',tnamstr,...
            'call','',...
            'visible','off',...
            'background',tcolor);
        ud.h.tval = uicontrol('style','edit',...
            'pos',[icvleft,tbot,icvalw,texth],...
            'string',tvalstr,...
            'horizon','left',...
            'user',1,...
            'call',solintcall,...
            'visible','off',...
            'background',ecolor);
        
        %% The figure
        
        setfight = eqfrbot + eqfrht +frsep;	% Height of the figure.
        setfigl = 30;
        setfigbot = 60;
        set(odesolveset,'pos',[setfigl setfigbot setfigw setfight]);
        
        
        
        ud.separation = separation;
        ud.left = left;
        ud.fudge = fudge;
        ud.frsep = frsep;
        ud.eqleft = eqleft;
        ud.varw = varw;
        ud.equalw = equalw;
        ud.tframeht = tframeht;
        ud.solfrbot = solfrbot;
        ud.icfrw = icfrw;
        ud.icfrl = icfrl;
        ud.frame = frame;
        ud.eqfrw = eqfrw;
        ud.solfrht = solfrht;
        ud.texth = texth;
        ud.eqlength = eqlength;
        ud.icnamew = icnamew;
        ud.icvalw = icvalw;
        ud.setfigw = setfigw;
        ud.ecolor = ecolor;
        ud.tcolor = tcolor;
        ud.dcolor = dcolor;
        ud.ccolor = ccolor;
        ud.dfname = '';
        ud.fileflag = 1;
        us.addflag = 0;
        ud.iflag = 1;    %% This is the initialization step.
        ud.settings.atol = .000001;
        ud.settings.rtol = .001;
        ud.settings.refine = 4;
        ud.settings.istep = 'computed';
        ud.settings.mstep = 'computed';
        ud.settings.ncval = 2;
        ud.defaultsettings = ud.settings;
        ud.solver = 'ode45';
        
        set(odesolveset,'user',ud);
        odesolve('numchange');
        ud = get(gcf,'user');
        NN = str2num(get(ud.h.NN,'string'));
        ud.NNN = 1:NN;
        for k = 1:NN
            set(ud.h.xname(k),'string',ud.c.xname{k});
            set(ud.h.der(k),'string',ud.c.der{k});
            set(ud.h.icname(k),'string',[ud.c.xname{k}, ' = ']);
            set(ud.h.icval(k),'string',ud.c.ic{k});
        end
        
        % The menus
        
        hhsetup = get(0,'showhiddenhandles');
        set(0,'showhiddenhandles','on');
        mefile = findobj(odesolveset,'label','&File');
        meedit = findobj(odesolveset,'label','&Edit');
        delete(findobj(odesolveset,'label','&Tools'));
        delete(findobj(odesolveset,'label','&View'));
        delete(findobj(odesolveset,'label','&Insert'));
        
        % The 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');
        meload = uimenu(mefile,'label','Load a system ...',...
            'call','odesolve(''loadsyst'',''system'');',...
            'pos',1);
        mesave = uimenu(mefile,'label','Save the current system ...',...
            'call','odesolve(''savesyst'',''system'');',...
            'pos',2);
        meloadg = uimenu(mefile,'label','Load a gallery ...',...
            'call','odesolve(''loadsyst'',''gallery'');',...
            'separator','on','pos',3);
        mesaveg = uimenu(mefile,'label','Save a gallery ...',...
            'call','odesolve(''savesyst'',''gallery'');',...
            'tag','savegal',...
            'pos',4);
        delgall = ['sud = get(gcf,''user'');',...
                'mh = get(sud.h.gallery,''children'');',...
                'add = findobj(sud.h.gallery,''tag'',''add system'');',...
                'mh(find(mh == add)) = [];',...
                'delete(mh);',...
                'set(sud.h.gallery,''user'',[]);',...
                'set(sud.h.loaddefgal,''enable'',''on'')'];
        
        
        medelg = uimenu(mefile,'label','Delete the current gallery.',...
            'call',delgall,'pos',5);
        melddg = uimenu(mefile,'label','Load the default gallery.',...
            'call','odesolve(''loadsyst'',''default'');',...
            'enable','off',...
            'tag','load default',...
            'pos',6);
        ud.h.loaddefgal = melddg;
        mesolnew = uimenu(mefile,...
            'label','Solve in a new window.',...
            'call',sncall,...
            'separator','on',...
            'accelerator','G','pos',7);
        mesoladd = uimenu(mefile,...
            'label','Add to the existing display.',...
            'call',socall,...
            'enable','off',...
            'separator','off',...
            'accelerator','H','pos',8);
        ud.h.somenu = mesoladd;
        set(mepset,'vis','on','pos',9);
        set(meprev,'vis','on','pos',10);
        set(meexp,'vis','on','pos',11,'separator','off');
        merestart = uimenu(mefile,'label',...
            'Restart ODESOLVE.',...
            'call','odesolve(''restart'')',...
            'separator','on','pos',12);
        
        mequit = uimenu(mefile,...
            'label','Quit ODESOLVE',...
            'call','odesolve(''quit'')',...
            'separator','off','pos',13);
        
        % The Edit menu
        
        set(get(meedit,'child'),'vis','off');
        eqclear = [
            'ud = get(gcf,''user'');h = ud.h;',...
                'set([h.xname,h.der,h.icname,h.icval],''string'','''');'];
        meclrf = uimenu(meedit,'label','Clear equations',...
            'call',eqclear,...
            'accelerator','E');
        pclear = [
            'ud = get(gcf,''user'');h = ud.h;',...
                'set([h.pname,h.pval],''string'','''');',...
                'ud.c.pname = {};',...
                'ud.c.pval = [];',...
                'set(gcf,''user'',ud);',...
            ];
        meclrp = uimenu(meedit,'label','Clear parameters',...
            'call',pclear,...
            'accelerator','P');
        exclear = [
            'ud = get(gcf,''user'');h = ud.h;',...
                'set([h.exname,h.exval],''string'','''');',...
                'ud.c.exname = {};',...
                'ud.c.exval = {};',...
                'set(gcf,''user'',ud);',...
            ];
        meclrex = uimenu(meedit,'label','Clear expressions',...
            'call',exclear,...
            'accelerator','X');
        allclear = [
            'ud = get(gcf,''user'');h = ud.h;',...
                'set([h.xname,h.der,h.icname,h.icval],''string'','''');',...
                'set([h.pname,h.pval],''string'','''');',...
                'set([h.exname,h.exval],''string'','''');',...
                'ud.c.pname = {};',...
                'ud.c.pval = [];',...
                'ud.c.exname = {};',...
                'ud.c.exval = {};',...
                'set(gcf,''user'',ud);',...
            ];
        meclrall = uimenu(meedit,'label','Clear all',...
            'call',allclear,...
            'accelerator','A',...
            'separator','on');
        
        % The Gallery Menu
        
        sysmenu = uimenu('label','Gallery','visible','off','pos',3);
        
        meadd = uimenu(sysmenu,'label','Add current system to the gallery',...
            'call','odesolve(''addgall'');',...
            'tag','add system');
        sep = 'on';
        for kk = 1:length(system)
            kkk = num2str(kk);
            if kk == 2, sep = 'off';end
            sysmen(kk) = uimenu(sysmenu,'label',system(kk).name,...
                'call',['odesolve(''system'',',kkk,')'],...
                'separator',sep,'visible','off');
        end
        set(sysmenu,'user',system);
        ud.h.gallery = sysmenu;
        ud.egg = (exist('EASTEREGG') ==2);
        
        % The Options Menu
        
        optmenu = uimenu('label','Options','visible','off','pos',4);
        ud.pplotflag = 0;
        pplotcall = [
            '[h,fig] = gcbo;',...
                'sud = get(fig,''user'');',...
                'sw = get(h,''checked'');',...
                'switch sw,',...
                ' case ''on'',',...
                '  set(h,''checked'',''off'');',...
                '  sud.pplotflag = 0;',...
                ' case ''off'',',...
                '  set(h,''checked'',''on'');',...
                '  sud.pplotflag = 1;',...
                'end,',...
                'set(fig,''user'',sud);',...
            ];
        
        
        
        mepplot = uimenu(optmenu,'label','Mark computed points.',...
            'call',pplotcall,...
            'visible','off');
        
        ud.icplotflag = 0;
        icplotcall = [
            '[h,fig] = gcbo;',...
                'sud = get(fig,''user'');',...
                'sw = get(h,''checked'');',...
                'switch sw,',...
                ' case ''on'',',...
                '  set(h,''checked'',''off'');',...
                '  sud.icplotflag = 0;',...
                ' case ''off'',',...
                '  set(h,''checked'',''on'');',...
                '  sud.icplotflag = 1;',...
                'end,',...
                'set(fig,''user'',sud);',...
            ];
        
        meicplot = uimenu(optmenu,'label','Mark initial points.',...
            'call',icplotcall,...
            'visible','off');
        
        mefileexport = uimenu(optmenu,'label','Export ODE function file.',...
            'call','odesolve(''fileexport'')',...
            'visible','off');
        
        mesolvset = uimenu(optmenu,'label','Solver settings.',...
            'call','odesolve(''settings'')',...
            'visible','off');
        
        
        
        % Record the handles in the User Data of the Set Up figure.
        
        set(odesolveset,'user',ud);
        hhhh = findobj(odesolveset,'type','uicontrol');
        set(hhhh,'units','normal')
        
        children = get(odesolveset,'children');
        ics = ud.h.icval;
        for kk = 1:NN
            children(find(children == ics(kk))) = [];
        end
        children = [fliplr(ics)';children];
        set(odesolveset,'children',children);
        
        set(odesolveset,'visible','on','resize','on');
        set(get(odesolveset,'children'),'visible','on');
        set(get(sysmenu,'children'),'visible','on');
        set(get(optmenu,'children'),'visible','on');
        set(0,'showhiddenhandles',hhsetup);
        
        
    case 'numchange'
        
        fig = findobj('name','ODESOLVE Setup');
        hhhh = findobj(fig,'type','uicontrol');
        set(hhhh,'units','pix')
        ud = get(fig,'user');
        iflag = ud.iflag;
        eqfrbot = ud.tframeht;
        separation = ud.separation;
        left = ud.left;
        fudge = ud.fudge;
        eqleft = ud.eqleft;
        varw = ud.varw;
        equalw = ud.equalw;
        tframeht = ud.tframeht;
        frame = ud.frame;
        eqfrw = ud.eqfrw;
        frsep = ud.frsep;
        texth = ud.texth;
        eqlength = ud.eqlength;
        icnamew = ud.icnamew;
        icvalw = ud.icvalw;
        ecolor = ud.ecolor;
        tcolor = ud.tcolor;
        dcolor = ud.dcolor;
        ccolor = ud.ccolor;
        xname = ud.h.xname;
        equals = ud.h.equals;
        der = ud.h.der;
        icname = ud.h.icname;
        icval = ud.h.icval;
        
        NN = str2num(get(ud.h.NN,'string'));  %% New number
        NNN = length(ud.h.xname);             %% Old number
        if NNN == NN        %% No changes needed.
            return
        elseif NNN > NN      %% Remove extra boxes and handles.
            delete(xname(NN+1:NNN));
            xname = xname(1:NN);
            ud.c.xname = ud.c.xname(1:NN);
            delete(equals(NN+1:NNN));
            equals = equals(1:NN);
            delete(der(NN+1:NNN));
            der = der(1:NN);
            ud.c.der = ud.c.der(1:NN);
            KK = length(ud.c.ic);
            if NN <= KK
                ud.c.ic = ud.c.ic(1:NN);
            end
            delete(icname(NN+1:NNN));
            icname = icname(1:NN);
            delete(icval(NN+1:NNN));
            icval = icval(1:NN);
        end
        
        eqfrht = (NN+1)*separation +10;
        eqfrwind = [left, eqfrbot, eqfrw, eqfrht];
        
        set(frame(1),'pos',eqfrwind);
        icfrbot = ud.solfrbot;
        icfrht = ud.solfrht + eqfrht +frsep;
        icfrw = ud.icfrw;
        icfrl = 2*left+eqfrw;
        icfrwind = [icfrl, icfrbot, icfrw, icfrht]; 
        set(frame(2),'pos',icfrwind);
        
        xnamecall = [ 
            '[h,fig] = gcbo;'...
                'ud =  get(fig,''user'');'...
                'num = get(h,''user'');'...
                'str = get(h,''string'');'...
                'ud.c.xname{num} = str;'...
                'set(ud.h.icname(num),''string'',[str, '' = '']);'...
                'ud.fileflag = 1;'...
                'set(fig,''user'',ud);' ];
        
        dercall = [
            '[h,fig] = gcbo;'...
                'ud =  get(fig,''user'');'...
                'num = get(h,''user'');'...
                'ud.c.der{num} = get(h,''string'');'...
                'ud.fileflag = 1;'...
                'set(h,''backgroundcolor'',ud.ecolor);'...
                'odesolve(''reset'');'...
                'set(fig,''user'',ud);' ];
        iccall = [
            '[h,fig] = gcbo;'...
                'ud =  get(fig,''user'');'...
                'value = get(h,''string'');'...
                'num = get(h,''user'');'...
                'ud.c.ic{num} = value;'...
                'set(h,''background'',ud.ecolor);'...
                'set(fig,''user'',ud);' ];
        
        equalsleft = ud.eqleft+varw;
        derleft = equalsleft+equalw;
        icnleft = icfrl +2;
        icvleft = icnleft + icnamew;
        for kk = 1:NN
            ebot(kk) = eqfrbot + 5 + (NN - kk)*separation;
        end   
        
        %% Move the boxes.
        
        KK = length(ud.c.xname);
        for kk = 1:min(NN,NNN)
            set(xname(kk),'pos',[eqleft, ebot(kk), varw, texth]);
            set(equals(kk),'pos',[equalsleft, ebot(kk)-fudge, equalw, texth]);
            set(der(kk),'pos',[derleft, ebot(kk), eqlength, texth]);
            set(icname(kk),'pos',[icnleft,ebot(kk)-fudge,icnamew,texth]);
            set(icval(kk),'pos',[icvleft,ebot(kk),icvalw,texth]);
        end
        
        %% Create new boxes if NN > NNN
        
        for kk =(NNN+1):NN
            if kk <= KK
                namstr = ud.c.xname{kk};
                derstr = ud.c.der{kk};
                icnamstr = [namstr, ' = '];
                icvalstr = ud.c.ic(kk); 
            else
                namstr = '';
                derstr = '';
                icnamstr = ' = ';
                icvalstr = ''; 
            end
            xname(kk) = uicontrol('style','edit',...
                'pos',[eqleft, ebot(kk), varw, texth],...
                'horizon','right',...
                'string',namstr,...
                'call',xnamecall,...
                'user',kk,...
                'visible','off',...
                'parent',fig,...
                'backgroundcolor',ecolor);
            ud.c.xname{kk} = namstr;
            equals(kk) = uicontrol('style','text',...
                'pos',[equalsleft,ebot(kk)-fudge,equalw,texth],...
                'horizon','center',...
                'string',''' = ',...
                'visible','off',...
                'parent',fig,...
                'backgroundcolor',tcolor);
            der(kk) = uicontrol('style','edit',...
                'pos',[derleft,ebot(kk),eqlength,texth],...
                'string',derstr,...
                'user',kk,...
                'horizon','left',...
                'call',dercall,...
                'visible','off',...
                'parent',fig,...
                'backgroundcolor',ecolor);
            ud.c.der{kk} = derstr;
            icname(kk) = uicontrol('style','text',...
                'pos',[icnleft,ebot(kk)-fudge,icnamew,texth],...
                'horizon','right',...
                'string',icnamstr,...
                'visible','off',...
                'parent',fig,...
                'background',tcolor);
            if iflag
                col = ecolor;
            else
                col = ecolor;
            end
            icval(kk) = uicontrol('style','edit',...
                'pos',[icvleft,ebot(kk),icvalw,texth],...
                'string',icvalstr,...
                'user',kk,...
                'call',iccall,...
                'horizon','left',...
                'visible','off',...
                'parent',fig,...
                'background',col);
        end
        eqlabbot = eqfrbot+5+NN*separation;
        eqlabw = varw+equalw+eqlength;
        set(ud.eqnlabel,'pos',[eqleft,eqlabbot,eqlabw,texth]);
        iclabbot = eqfrbot+5+NN*separation;
        iclabl = ud.icfrl+left;
        iclabw = icfrw - 2*left;
        set(ud.iclabel,'pos',[iclabl,iclabbot,iclabw,texth]);
        ud.c.NN = NN;
        ud.h.xname = xname;
        ud.h.equals = equals;
        ud.h.der = der;
        ud.h.icname = icname;
        ud.h.icval = icval;
        ud.iflag = 0;
        ud.NNN = (1:length(xname))';
        setfight = eqfrbot + eqfrht + frsep;	% Height of the figure.
        setfigl = 30;
        setfigbot = 60;
        set(fig,'pos',[setfigl setfigbot ud.setfigw setfight])
        set(fig,'vis','on','user',ud);
        odesolve('reset');
        set(get(fig,'children'),'vis','on');  
        
        hhhh = findobj(fig,'type','uicontrol');
        set(hhhh,'units','normal')
        if NN == 1
            set(ud.h.output,'enable','off')
        else
            set(ud.h.output,'enable','on');
        end
        
    case 'reset'
        
        odset = findobj('name','ODESOLVE Setup');
        ud = get(odset,'user');
        set(ud.h.output,'value',1)
        str = ['     all variables vs ', ud.c.tname, '.'];
        set(ud.h.outputstr,'string',str);
        ud.outputhandle = @tpout;
        set(ud.h.sobutt,'enable','off');
        if isfield(ud.h,'somenu')
            set(ud.h.somenu,'enable','off');
        end
        set(odset,'user',ud);
        
    case 'savesyst'
        
        odset = findobj('name','ODESOLVE Setup');
        type = input1;
        sud = get(odset,'user');
        systems = get(sud.h.gallery,'user');
        
        switch type
            case 'system'  
                fn = sud.h.name;
                if ~isempty(fn)
                    fn(find(abs(fn)==32))='_';   % Replace spaces by underlines.
                end	
                fname = [fn,'.ods'];
                comp = computer;
                switch  comp
                    case 'PCWIN'
                        filter = [sud.oddir,'\',fname];
                    case 'MAC2'
                        filter = [sud.oddir,':',fname];
                    otherwise
                        filter = [sud.oddir,'/',fname];
                end
                [fname,pname] = uiputfile(filter,'Save the system as:');
                if fname == 0,return;end
                ll = length(fname);
                if (ll>4 & strcmp(fname(ll-3:ll),'.ods'))
                    fn = fname(1:(ll-4));
                else
                    fn =fname;
                end
                fname = [fn,'.ods'];
                sud.h.name = fn;
                set(odset,'user',sud);
                sud.c.name = fn;
                newsysts = sud.c;
                
            case 'gallery'
                ll = length(systems); 
                if ll == 0
                    warndlg(['There are no systems to make up a gallery.'],'Warning');
                    return
                end
                [names{1:ll}] = deal(systems.name);
                [sel,ok] = listdlg('PromptString','Select the systems',...
                    'Name','Gallery selection',...
                    'ListString',names);
                if isempty(sel)
                    return
                else
                    newsysts = systems(sel);
                end	
                comp = computer;
                switch  comp
                    case 'PCWIN'
                        prompt = [sud.oddir,'\*.odg'];
                    case 'MAC2'
                        prompt = [sud.oddir,':*.odg'];
                    otherwise
                        prompt = [sud.oddir,':/.odg'];
                end
                [fname,pname] = uiputfile(prompt,'Save the gallery as:');
                if fname == 0,return;end
                ll = length(fname);
                if (ll>4 & strcmp(fname(ll-3:ll),'.odg'))
                    fn = fname(1:(ll-4));
                else
                    fn =fname;
                end
                fname = [fn,'.odg'];
            end
            
            ll = length(newsysts);
            fid = fopen([pname fname],'w');
            for k = 1:ll
                NN = newsysts(k).NN;
                nstr = ['H.name = ''', newsysts(k).name, ''';\n'];
                fprintf(fid,nstr);
                NNstr = ['H.NN = ', num2str(NN), ';\n'];
                fprintf(fid,NNstr);
                xname = newsysts(k).xname;
                xnstr = ['H.xname = {''', xname{1}, ''''];
                for kk = 2:NN
                    xnstr = [xnstr, ',''', xname{kk}, ''''];
                end
                xnstr = strrep(xnstr,'\','\\');
                xnstr = [xnstr, '};\n'];
                fprintf(fid,xnstr);
                der = strrep(newsysts(k).der,'\','\\');
                derstr = ['H.der = {''', der{1}, ''''];
                for kk = 2:NN
                    derstr = [derstr, ',''', der{kk}, ''''];
                end
                derstr = [derstr, '};\n'];
                fprintf(fid,derstr);
                ics = newsysts(k).ic;
                icstr = ['H.ic = {''', ics{1}, ''''];
                for kk = 2:NN
                    icstr = [icstr, ',''',ics{kk}, ''''];
                end
                icstr = [icstr, '};\n'];
                fprintf(fid,icstr);
                tnstr = ['H.tname = ''', newsysts(k).tname, ''';\n'];
                fprintf(fid,tnstr);
                tvstr =['H.tvals = [', num2str(newsysts(k).tvals), '];\n'];
                fprintf(fid,tvstr);
                NNN = length(newsysts(k).pname);
                if NNN == 0
                    pnstr = 'H.pname = {};\n';
                    pvstr = 'H.pval = [];\n';
                else
                    pname = strrep(newsysts(k).pname,'\','\\');
                    pnstr = ['H.pname = {''', pname{1}, ''''];
                    for kk = 2:NNN
                        pnstr = [pnstr, ',''', pname{kk}, ''''];
                    end
                    pnstr = [pnstr, '};\n'];
                    pvstr = ['H.pval = [', num2str(newsysts(k).pval), '];\n'];
                end
                fprintf(fid,pnstr);
                fprintf(fid,pvstr);
                NNN = length(newsysts(k).exname);
                if NNN == 0
                    exnstr = 'H.exname = {};\n';
                else
                    exname = strrep(newsysts(k).exname,'\','\\');
                    exnstr = ['H.exname = {''', exname{1}, ''''];
                    for kk = 2:NNN
                        exnstr = [exnstr, ',''', exname{kk}, ''''];
                    end
                    exnstr = [exnstr, '};\n'];
                end
                fprintf(fid,exnstr);
                NNN = length(newsysts(k).exval);
                if NNN == 0
                    exvstr = 'H.exval = {};\n';
                else
                    exval = strrep(newsysts(k).exval,'\','\\');
                    exvstr = ['H.exval = {''', exval{1}, ''''];
                    for kk = 2:NNN
                        exvstr = [exvstr, ',''', exval{kk}, ''''];
                    end
                    exvstr = [exvstr, '};\n'];
                end
                fprintf(fid,exvstr);
            end
            fclose(fid);
            
        case 'loadsyst'
            
            sud = get(gcf,'user');
            pos = get(gcf,'pos');
            wpos = [pos(1),pos(2)+pos(4)+20,300,20];
            waith = figure('pos',wpos,...
                'numb','off',...
                'vis','off',...
                'next','replace',...
                'menubar','none',...
                'resize','off',...
                'createfcn','');
            axes('pos',[0.01,0.01,0.98,0.98],...
                'vis','off');
            xp = [0 0 0 0];
            yp = [0 0 1 1];
            xl = [1 0 0 1 1];
            yl = [0 0 1 1 0];
            patchh = patch(xp,yp,'r','edgecolor','r','erase','none');
            lineh = line(xl,yl,'erase','none','color','k');
            type = input1;  
            set(sud.h.gallery,'enable','off');
            
            switch type
                case 'default'
                    set(waith,'name','Loading the default gallery.','vis','on');
                    set(findobj('tag','load default'),'enable','off');
                    megall = sud.h.gallery;
                    mh = get(megall,'children');
                    add = findobj(megall,'tag','add system');
                    mh(find(mh == add)) = [];
                    delete(mh);
                    newsysstruct = get(megall,'user');
                    system = sud.system;
                    LL = length(system);
                    x = 1/(LL+2);
                    xp = [xp(2),x,x,xp(2)];
                    set(patchh,'xdata',xp);
                    set(lineh,'xdata',xl);
                    drawnow;
                    sep = 'on';
                    for kk = 1:length(system)
                        kkk = num2str(kk);
                        if kk == 2, sep = 'off';end
                        uimenu(megall,'label',system(kk).name,...
                            'call',['odesolve(''system'',',kkk,')'],...
                            'separator',sep);
                    end % for
                    set(megall,'user',system);
                    set(sud.h.loaddefgal,'enable','off')
                    
                otherwise
                    comp = computer;
                    switch  comp
                        case 'PCWIN'
                            prompt = [sud.oddir,'\'];
                        case 'MAC2'
                            prompt = [sud.oddir,':'];
                        otherwise
                            prompt = [sud.oddir,'/'];
                    end
                    
                    switch type
                        case 'system'
                            prompt = [prompt,'*.ods'];
                            [fname,pname] = uigetfile(prompt,'Select a system to load.');
                        case 'gallery'
                            prompt = [prompt,'*.odg'];
                            [fname,pname] = uigetfile(prompt,'Select a gallery to load.');
                    end  % switch type
                    if fname == 0
                        delete(waith);
                        set(sud.h.gallery,'enable','on');
                        return;
                    end
                    set(waith,'name',['Loading ',fname],'vis','on');
                    fid = fopen([pname fname],'r');
                    newsysts = {};
                    kk = 0;
                    while ~feof(fid)
                        kk = kk + 1;
                        newsysts{kk} = fgetl(fid);
                    end
                    fclose(fid);
                    if mod(kk,11)
                        switch type
                            case 'system'
                                warndlg(['The file ',fname, ' does not define a proper system.'],...
                                    'Warning');
                            case 'gallery'
                                warndlg(['The file ',fname, ' does not define a proper gallery.'],...
                                    'Warning');
                        end
                        set(sud.h.gallery,'enable','on');
                        delete(waith);
                        return
                    end % if mod(kk,11)
                    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 = 1:11;
                            eval(newsysts{(j-1)*11+k});
                        end
                        newsysstruct(j) = H;
                    end
                    %    end
                    ignoresyst = {};
                    for j = 1:nnn
                        x = (j+1)/(nnn+2);
                        xp = [xp(2),x,x,xp(2)];
                        set(patchh,'xdata',xp);
                        set(lineh,'xdata',xl);
                        drawnow;
                        newsyst = newsysstruct(j);
                        sname = newsyst.name;
                        sname(find(abs(sname) == 95)) = ' '; % Replace underscores with spaces.
                        newsyst.name = sname;
                        ignore = odesolve('addgall',newsyst);
                        if ignore == -1;
                            ignoresyst{length(ignoresyst)+1} = sname;
                        end  
                    end
                    l = length(ignoresyst);
                    if l  % There was at least one system which was a dup with a different name.
                        if l == 1
                            message = {['The system "',ignoresyst{1},'" duplicates a ',...
                                        'system already in the gallery and was not added.']};
                        else
                            message = 'The systems ';
                            for k = 1:(l-1)
                                message = [message,'"',ignoresyst{k},'", ']; 
                            end
                            message = {[message,'and "',ignoresyst{l},'" duplicate ',...
                                        'systems already in the gallery and were not added.']};
                        end % if l == 1 & else
                        helpdlg(message,'Ignored systems');
                    end  % if l
                    if strcmp(type,'system') % Added a system.
                        if ignore > 0 % The system was ignored.
                            kk = ignore;
                        else
                            systems = get(sud.h.gallery,'user');
                            kk = length(systems);
                        end
                        odesolve('system',kk);
                    end 
            end
            if strcmp('type','default')
                odesolve('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);
            
            
        case 'addgall'
            
            output = 0;
            odset = findobj('name','ODESOLVE Setup');
            sud = get(odset,'user');
            if nargin < 2    % We are adding the current system.
                syst = sud.c;
                sname = inputdlg('Provide a name for this system.','System name',1,{syst.name});
                if isempty(sname),return;end
                sname = sname{1};
                if ~strcmp(sname,syst.name)
                    sud.c.name = sname;
                    set(odset,'user',sud);
                    syst.name = sname;
                end  
            else  % We have a system coming from a file.
                syst = input1;
                sname = syst.name;
            end
            systems = get(sud.h.gallery,'user');
            LL = length(systems);
            kk = 1;
            while ((kk<=LL) & (~strcmp(sname,systems(kk).name)))
                kk = kk + 1;
            end
            nameflag = (kk<=LL);
            ssyst = rmfield(syst,'name');
            kk = 1;
            while ((kk<=LL) & (~isequal(ssyst,rmfield(systems(kk),'name'))))
                kk = kk + 1;
            end
            systflag = 2*(kk<=LL);
            flag = nameflag + systflag;
            switch flag
                case 1  % Same name but different system.
                    mh = findobj(sud.h.gallery,'label',sname);
                    prompt = {['The system "',sname,'", which you wish to add to the gallery has ',...
                                'the same name as a different system already in the gallery.  Please ',...
                                'specify the name for the newly added system.'],...
                            'Specify the name for the old system.'};
                    title = 'Two systems with the same name';
                    lineno = 1;
                    defans = {sname,sname};
                    answer = inputdlg(prompt,title,lineno,defans);
                    if isempty(answer)
                        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,'user',systems);
                syst.name = sname;
                if flag <=1
                    switch LL
                        case 0
                            systems = syst;
                            sepstr = 'on';
                        case 8
                            systems(9) = syst;
                            if strcmp(systems(8).name,'Duffing''s equation')
                                sepstr = 'on';
                            else
                                sepstr = 'off';
                            end
                        otherwise
                            systems(LL+1) = syst;
                            sepstr = 'off';
                    end      
                    kkk = num2str(LL+1);
                    newmenu = uimenu(sud.h.gallery,'label',sname,...
                        'call',['odesolve(''system'',',kkk,')'],...
                        'separator',sepstr);
                    set(findobj('tag','savegal'),'enable','on');
                end
                set(sud.h.gallery,'user',systems);
                
            case 'system'
                
                odset = findobj('name','ODESOLVE Setup');
                ud = get(odset,'user');
                ud.fileflag = 1;
                kk = input1;
                if isstr(kk)
                    kk = str2num(input1);
                end
                system = get(ud.h.gallery,'user');
                syst = system(kk);
                ud.o = syst;
                ud.c = syst;
                ud.h.name = syst.name;
                NN = syst.NN;
                str = ['     all variables vs ', syst.tname, '.'];
                set(ud.h.outputstr,'string', str);
                ud.NNN = 1:NN;
                set(ud.h.NN,'string',num2str(syst.NN));
                set(odset,'user',ud);
                odesolve('numchange');
                ud = get(odset,'user');
                for k = 1:NN
                    set(ud.h.xname(k),'string',syst.xname{k},'vis','on');
                    set(ud.h.der(k),'string',syst.der{k});
                    set(ud.h.icname(k),'string',[syst.xname{k}, ' = ']);
                    set(ud.h.icval(k),'string',syst.ic{k});
                end
                tname = syst.tname;
                pname = syst.pname;
                pval = syst.pval;
                
                for k = 1:length(ud.h.pname)
                    if k > length(pname)
                        pname{k} = '';
                    end
                    name = pname{k};
                    if k <= length(pval)
                        value = num2str(pval(k));
                    else
                        value = '';
                    end
                    set(ud.h.pname(k),'string',name);
                    set(ud.h.pval(k),'string',value);
                end
                
                exval = syst.exval;
                exname = syst.exname;
                for k = 1:length(ud.h.exname)
                    if k > length(exname)
                        exname{k} = '';
                    end
                    name = exname{k};
                    if k <= length(exval)
                        value = num2str(exval{k});
                    else
                        value = '';
                    end
                    set(ud.h.exname(k),'string',name);
                    set(ud.h.exval(k),'string',value);
                end
                ud.c.exname = exname;
                ud.c.exval = exval;
                ud.c.pname = pname;
                ud.c.pval = pval;
                tn = tname;
                set(ud.h.tvar,'string',tn);
                set(ud.h.solint(1),'string',[' <=  ', tn, '  <= ']);
                str = num2str(syst.tvals(2));
                set(ud.h.solint(2),'string',str);
                str = num2str(syst.tvals(3));
                set(ud.h.solint(3),'string',str);
                set(ud.h.tic,'string',[tn, ' = ']);
                set(ud.h.tval,'string',num2str(syst.tvals(1)));
                set(odset,'user',ud);
                odesolve('reset');
                
            case 'solve'
                
                %%  Change the defaults.  We might want to delay this a little.
                
                odesolveset = findobj('name','ODESOLVE Setup');
                sud = get(odesolveset,'user');
                sud.o = sud.c;
                set(odesolveset,'user',sud);
                
                %% Get the pertinent information
                
                NN = sud.c.NN;
                xname = sud.c.xname;
                tname = sud.c.tname;
                der = sud.c.der;
                % ic = sud.c.ic;
                tvals = sud.c.tvals;
                tinit = tvals(1);
                solvers = get(sud.h.solver,'string');
                value = get(sud.h.solver,'value');
                solver = str2func(solvers{value});
                
                %% Check that all of the ICs have an entry.
                
                hempt = findobj(sud.h.icval,'string','');
                if ~isempty(hempt)
                    set(hempt,'backgroundcolor',sud.dcolor);
                    mtitle = 'Initial condition error';
                    mmsg = 'Some initial conditions have no values.';
                    msgbox(mmsg,mtitle);
                    return
                end
                
                % Find the parameters and expressions that are actually used.
                
                [renn,rpnn] = pandecheck('write');
                
                exname = sud.h.exname;
                exval = sud.h.exval;
                krex = length(renn);
                rexname = cell(1,krex);
                for kk = 1:krex
                    rexname{kk} = get(exname(renn(kk)),'string');
                    rexval{kk} = get(exval(renn(kk)),'string');
                end
                
                pname = sud.h.pname;
                pval = sud.h.pval;
                icond = sud.c.ic;
                krp = length(rpnn);
                rpval = cell(1,krp);
                for kk = 1:krp
                    pn = get(pname(rpnn(kk)),'string');
                    rpname{kk} = [',', pn];
                    rpns{kk} = pn;
                    rpval{kk} = get(pval(rpnn(kk)),'string');
                end  
                
                %% Check that the parameter values are numbers.
                
                if ~isempty(rpnn)
                    hbad = [];
                    for kk = 1:length(rpnn)
                        if isempty(str2num(rpval{kk}))
                            hbad = [hbad, pval(rpnn(kk))];
                        end
                    end
                    if ~isempty(hbad)
                        set(hbad ,'backgroundcolor',sud.dcolor);
                        mtitle = 'Parameter error';
                        mmsg = 'Some parameters do not have numerical values.';
                        msgbox(mmsg,mtitle,'modal');
                        return
                    end      
                end
                
                rpvalstr = cell(1,0);
                for k = 1:length(rpval)
                    rpvalstr{k} = str2num(rpval{k});
                end
                
                
                %% Compute the initial conditions, which may involve parameters.
                
                ic = zeros(NN,1);
                hbad = [];
                for kk = 1:NN           % Replace parameters by their values.
                    str = deblank(sud.c.ic{kk});
                    if exist('rpname')
                        lrp =length(rpns);
                        for jj = 1:lrp
                            str = pareval(str,rpns{jj},rpval{jj});
                        end
                    end
                    try
                        ic(kk) = eval(str);
                    catch
                        hbad = [hbad,sud.h.icval(kk)];
                    end
                end
                if ~isempty(hbad )
                    set(hbad ,'backgroundcolor',sud.dcolor);
                    mtitle = 'Initial condition error';
                    mmsg = 'Some initial conditions do not evaluate to numerical values.';
                    msgbox(mmsg,mtitle,'modal');
                    return
                end
                
                %% Compute the solution interval.  The endpoints may involve
                %%  parameters. 
                
                % Yet to be done.
                
                
                % Error check the expressions.
                
                if krex
                    hbad = [];
                    for kk = 1:krex
                        str = rexval{kk};
                        str = pareval(str,tname,num2str(tinit));
                        for jj = 1:krp
                            str = pareval(str,rpns{jj},rpval{jj});
                        end
                        for jj = 1:NN
                            str = pareval(str,xname{jj},num2str(ic(jj)));
                        end
                        try
                            exvalue(kk) = eval(str);
                        catch
                            hbad = [hbad,exval(renn(kk))];
                        end
                    end
                    if ~isempty(hbad)
                        set(hbad ,'backgroundcolor',sud.dcolor);
                        mtitle = 'Expression error';
                        mmsg = 'Some expressions are not entered correctly.';
                        msgbox(mmsg,mtitle,'modal');
                        return
                    end      
                end
                
                % Error check the derivative strings.
                
                hbad = [];
                for kk = 1:NN
                    str = der{kk};
                    str = pareval(str,tname,num2str(tinit));
                    for jj = 1:krp
                        str = pareval(str,rpns{jj},rpval{jj});
                    end
                    for jj = 1:NN
                        str = pareval(str,xname{jj},num2str(ic(jj)));
                    end
                    for jj = 1:krex
                        str = pareval(str,rexname{jj},num2str(exval(jj)));
                    end
                    try
                        derval(kk) = eval(str);
                    catch
                        hbad = [hbad, sud.h.der(kk)];
                    end
                end
                if ~isempty(hbad)
                    set(hbad ,'backgroundcolor',sud.dcolor);
                    mtitle = 'Equation error';
                    mmsg = 'Some equations are not entered correctly.';
                    msgbox(mmsg,mtitle,'modal');
                    return
                end      
                
                % Start a new derivative m-file if necessary.
                
                dfcn = sud.dfname;
                if sud.fileflag
                    if exist(dfcn) == 2
                        delete([tempdir, dfcn, '.m']);
                    end
                    
                    tee = clock;
                    tee = ceil(tee(6)*100);
                    dfcn=['odestp',num2str(tee)];
                    
                    writefile(dfcn, tempdir, renn, rpnn);
                    
                    sud.dfname = dfcn;
                    sud.fileflag = 0;
                    set(odesolveset,'user',sud);
                end
                
                %% Solve the system.
                
                intplus = tvals([1 3]);
                intminus = tvals([1 2]);
                settings = sud.settings;
                opt = odeset('outputfcn',sud.outputhandle,...
                    'abstol',settings.atol,...
                    'reltol',settings.rtol,...
                    'refine',settings.refine);
                if ~isstr(settings.istep)
                    opt = odeset(opt,'initialstep',settings.istep);
                end
                if ~isstr(settings.mstep)
                    opt = odeset(opt,'maxstep',settings.mstep);
                end
                if settings.ncval == 1
                    opt = odeset(opt,'normcontrol','on');
                end
                
                odesdisp = findobj('name','ODESOLVE Display');
                NNN = sud.NNN;
                
                if isempty(odesdisp)
                    odesdisp = figure('name','ODESOLVE Display',...
                        'numb','off',...
                        'tag','odesolve');
                    odesolve('figdefault',odesdisp);
                    dud.solutions = [];
                    dud.axes = axes('units','pix',...
                        'box','on',...
                        'interrupt','on',...
                        'xgrid','on',...
                        'ygrid','on',...
                        'drawmode','fast',...
                        'tag','display axes');
                    set(sud.h.sobutt,'enable','on');
                    set(sud.h.somenu,'enable','on');
                    
                    doptmenu = uimenu('label','Options','visible','off','pos',3);
                    
                    meexportdata = uimenu(doptmenu,'label','Export solution data.',...
                        'call','odesolve(''export'')',...
                        'separator','off','tag','dexp');
                    
                    
                    sud.addflag = 0;
                    sud.h.disp = odesdisp;
                    set(odesolveset,'user',sud);
                    set(odesdisp,'user',dud);
                    set(get(odesdisp,'children'),'vis','on')
                    set(get(doptmenu,'children'),'vis','on')
                else
                    figure(sud.h.disp);
                end
                
                dud = get(odesdisp,'user');
                dud.ic = ic;
                set(odesdisp,'user',dud);
                if sud.addflag
                    set(dud.axes,'next','add');
                else
                    set(dud.axes,'next','replace');
                    dud.solutions = [];
                    set(odesdisp,'user',dud);
                end
                
                exist(dfcn);
                dfh = str2func(dfcn);
                cflag = 0;
                
                if intplus(1) < intplus(2)
                    cflag = cflag + 1;
                    [tp,xp] = feval(solver,dfh,intplus,ic,opt,rpvalstr{:});
                    dud = get(odesdisp,'user');
                    hnew1 = dud.lines;
                    set(dud.axes,'next','add');
                end
                if intminus(1) > intminus(2)
                    cflag = cflag + 2;
                    [tm,xm] = feval(solver,dfh,intminus,ic,opt,rpvalstr{:});
                    dud = get(odesdisp,'user');
                    hnew2 = dud.lines;
                end
                
                dud.lines = [];
                set(odesdisp,'user',dud);
                
                % Store the trajectory.
                
                switch cflag
                    case 1 % positive only
                        x = xp;
                        t = tp;
                        dud.currentplot = hnew1;
                        
                    case 2 % negative only
                        x = flipud(xm);
                        t = flipud(tm);
                        dud.currentplot = hnew2;
                        
                    case 3 % both directions
                        x = flipud(xm);
                        t = flipud(tm);
                        x=[x;xp];
                        t=[t;tp];
                        delete(hnew2)
                        dud.currentplot = hnew1;
                        
                end	 % switch cflag
                
                LL = length(dud.solutions);
                if LL == 0
                    dud.solutions = struct('tdata',t,'xdata',x,'handles',dud.currentplot);
                else
                    dud.solutions(LL+1).tdata = t;
                    dud.solutions(LL+1).xdata = x;
                    dud.solutions(LL+1).handles = dud.currentplot;
                end
                set(dud.axes,'next','replace');
                set(odesdisp,'user',dud);
                set(dud.axes,'user',LL+1);
                val = get(sud.h.output,'value');
                set(sud.h.sobutt,'enable','on');
                set(sud.h.somenu,'enable','on');
                switch val
                    case 1
                        odesolve('tpdone');
                    case 2
                        odesolve('ppdone');
                    case 3
                        odesolve('p3done');
                end % switch val
                
            case 'tpdone'
                
                odesolveset = findobj('name','ODESOLVE Setup');
                sud = get(odesolveset,'user');
                odesdisp = findobj('name','ODESOLVE Display');
                dud = get(odesdisp,'user');
                NNN = sud.NNN;
                tvals = sud.c.tvals;
                axh = dud.axes;
                solution = dud.solutions(get(axh,'user'));
                t = solution.tdata;
                x = solution.xdata;
                
                
                x = x(:,NNN);
                N = length(t);
                if (tvals(3) == t(N) ) & (tvals(2) == t(1))
                    uu = sort(x(:));
                    NN = length(uu);
                    del = (uu(NN) - uu(1))/10;
                    ymin = uu(1)-del;
                    ymax = uu(NN)+del;
                else 
                    tf = tvals(1) +(t(N) - tvals(1))/2;
                    ti = tvals(1) +(t(1) - tvals(1))/2;
                    K = find( (t<tf)&(t>ti) );
                    uu = x(K,:);
                    uu = sort(uu(:));
                    NN = length(uu);
                    rr = 0.1;
                    n1 = ceil(rr*NN);
                    n2 = floor((1-rr)*NN);
                    del = uu(n2) - uu(n1);
                    if del == 0
                        del = 1;
                    end
                    ymin = max(uu(1) - rr*del, uu(n1) - del);
                    ymin = min(ymin,min(dud.ic));
                    ymax = min(uu(NN) + rr*del, uu(n2) +del);
                    ymax = max(ymax,max(dud.ic));
                end
                if sud.addflag
                    ylim = get(dud.axes,'ylim');
                    ymin = min(ymin,ylim(1));
                    ymax = max(ymax,ylim(2));
                end
                for j = 1:length(NNN)
                    set(dud.currentplot(j),'xdata',t,'ydata',x(:,j));
                end
                if sud.icplotflag
                    hold on
                    plot(tvals(1),dud.ic(NNN),'k.','markersize',6);
                    hold off
                end
                set(dud.axes,'ylim',[ymin,ymax]);
                grid on
                xlabel(sud.c.tname)
                xname = sud.c.xname;
                xname = xname(NNN);
                LL = length(xname);
                ystr = xname{1};
                if LL == 2
                    ystr = [ystr, ' and ' xname{2}];
                elseif LL > 2
                    for j = 2:[LL-1]
                        ystr = [ystr, ', ', xname{j}];
                    end
                    ystr = [ystr, ', and ' xname{LL}];
                end
                ylabel(ystr)
                legend(xname,0)
                
            case 'ppdone'
                
                odesolveset = findobj('name','ODESOLVE Setup');
                sud = get(odesolveset,'user');
                odesdisp = findobj('name','ODESOLVE Display');
                dud = get(odesdisp,'user');
                NNN = sud.NNN;
                tvals = sud.c.tvals;
                axh = dud.axes;
                solution = dud.solutions(get(axh,'user'));
                t = solution.tdata;
                x = solution.xdata;
                x = x(:,NNN);
                N = length(t);
                if (tvals(3) == t(N) ) & (tvals(2) == t(1))  
                    
                    % If the solution is over the complet t interval choose the figure
                    % limits to show all of the solution and more.
                    uu = sort(x(:,1));
                    del = (uu(N) - uu(1))/10;
                    if del == 0
                        del = 1;
                    end
                    xmin = uu(1)-del;
                    xmax = uu(N)+del;
                    vv = sort(x(:,2));
                    del = (vv(N) - vv(1))/10;
                    if del == 0
                        del = 1;
                    end
                    ymin = vv(1)-del;
                    ymax = vv(N)+del;
                else 
                    
                    % Otherwise limit ourselves to 1/2 of the limits in both
                    % directions. However, include the initial point for sure. 
                    
                    tf = tvals(1) +(t(N) - tvals(1))/2;
                    ti = tvals(1) +(t(1) - tvals(1))/2;
                    K = find( (t<tf)&(t>ti) );
                    uu = x(K,1); vv = x(K,2);
                    uu = sort(uu);  vv = sort(vv);
                    NN = length(uu);
                    rr = 0.1;
                    n1 = ceil(rr*NN);
                    n2 = floor((1-rr)*NN);
                    ic = dud.ic(NNN);
                    del = uu(n2) - uu(n1);
                    if del == 0
                        del = 1;
                    end
                    xmin = max(uu(1) - rr*del, uu(n1) - del);
                    xmin = min(xmin,ic(1));
                    xmax = min(uu(NN) + rr*del, uu(n2) + del);
                    xmax = max(xmax,ic(1));
                    del = vv(n2) - vv(n1);
                    if del == 0
                        del = 1;
                    end
                    ymin = max(vv(1) - rr*del, vv(n1) - del);
                    ymin = min(ymin,ic(2));
                    ymax = min(vv(NN) + rr*del, vv(n2) + del);
                    ymax = max(ymax,ic(2));
                end
                if sud.addflag
                    xlim = get(dud.axes,'xlim');
                    xmin = min(xmin,xlim(1));
                    xmax = max(xmax,xlim(2));
                    ylim = get(dud.axes,'ylim');
                    ymin = min(ymin,ylim(1));
                    ymax = max(ymax,ylim(2));
                end
                set(dud.currentplot,'xdata',x(:,1),'ydata',x(:,2));
                if sud.icplotflag
                    hold on
                    plot(dud.ic(NNN(1)),dud.ic(NNN(2)),'k.','markersize',6)
                    hold off
                end
                set(dud.axes,'xlim',[xmin,xmax],'ylim',[ymin,ymax]);
                grid on
                xlabel(sud.c.xname(NNN(1)));
                ylabel(sud.c.xname(NNN(2)));
                
            case 'p3done'
                
                odesolveset = findobj('name','ODESOLVE Setup');
                sud = get(odesolveset,'user');
                odesdisp = findobj('name','ODESOLVE Display');
                dud = get(odesdisp,'user');
                NNN = sud.NNN;
                NN = sud.c.NN;
                tvals = sud.c.tvals;
                axh = dud.axes;
                solution = dud.solutions(get(axh,'user'));
                t = solution.tdata;
                x = solution.xdata;
                kt = 0;
                for k = 1:3
                    if NNN(k) > NN
                        X(:,k) = t;
                        Y(k) = tvals(1);
                        kt = k;  % This is the t-variable.
                    else
                        X(:,k) = x(:,NNN(k));
                        Y(k) = dud.ic(NNN(k));
                    end
                end
                N = length(t);
                if (tvals(3) == t(N) ) & (tvals(2) == t(1))  
                    
                    % If the solution is over the complet t interval choose the figure
                    % limits to show all of the solution and more.
                    xmin = zeros(1,3);
                    xmax = zeros(1,3);
                    for k = 1:3
                        uu = sort(X(:,k));
                        del = (uu(N) - uu(1))/10;
                        xmin(k) = uu(1)-del;
                        xmax(k) = uu(N)+del;
                    end
                else 
                    
                    
                    % Otherwise limit ourselves to 1/2 of the limits in both
                    % directions. However, include the initial point for sure. 
                    
                    tf = tvals(1) +(t(N) - tvals(1))/2;
                    ti = tvals(1) +(t(1) - tvals(1))/2;
                    K = find( (t<tf)&(t>ti) );
                    NK = length(K);
                    rr = 0.1;
                    n1 = ceil(rr*NK);
                    n2 = floor((1-rr)*NK);
                    for k = 1:3
                        if k == kt
                            xmin(k) = t(1);
                            xmax(k) = t(N); 
                        else 
                            uu = X(K,k);
                            uu = sort(uu);
                            del = uu(n2) - uu(n1);
                            xmin(k) = max(uu(1) - rr*del, uu(n1) - del);
                            ic = dud.ic(NNN(k));
                            xmin(k) = min(xmin(k),ic);
                            xmax(k) = min(uu(NN) + rr*del, uu(n2) +del);
                            xmax(k) = max(xmax(k),ic);
                        end
                    end
                end
                if sud.addflag
                    xlim = get(dud.axes,'xlim');
                    xmin(1) = min(xmin(1),xlim(1));
                    xmax(1) = max(xmax(1),xlim(2));
                    ylim = get(dud.axes,'ylim');
                    xmin(2) = min(xmin(2),ylim(1));
                    xmax(2) = max(xmax(2),ylim(2));
                    zlim = get(dud.axes,'zlim');
                    xmin(3) = min(xmin(3),zlim(1));
                    xmax(3) = max(xmax(3),zlim(2));
                end
                for k = 1:3
                    if xmin(k) == xmax(k)
                        if xmin(k) == 0
                            del = 1;
                        else
                            del = xmin(k)/10;
                        end
                        xmin(k) = xmin(k)-del;
                        xmax(k) = xmax(k)+del;
                    end
                end
                set(dud.currentplot,'xdata',X(:,1),...
                    'ydata',X(:,2),...
                    'zdata',X(:,3));
                set(dud.axes,'xlim',[xmin(1),xmax(1)],...
                    'ylim',[xmin(2),xmax(2)],...
                    'zlim',[xmin(3),xmax(3)]);
                if sud.icplotflag
                    hold on
                    plot3(Y(1),Y(2),Y(3),'k.','markersize',6)
                    hold off
                end
                grid on
                for k = 1:3
                    if k == kt
                        name{k} = sud.c.tname;
                    else
                        name{k} = sud.c.xname{NNN(k)};
                    end
                end
                xlabel(name(1));
                ylabel(name(2));
                zlabel(name(3));
                
                
            case 'figdefault'
                
                fig = input1;
                set(fig,'CloseRequestFcn','odesolve(''closefcn'')');
                dfset = findobj('name','ODESOLVE Setup');
                sud = get(dfset,'user');
                ud = get(fig,'user');
                ud.ssize = sud.ssize;
                fs = sud.fontsize;
                ud.fontsize = fs;
                set(fig,'defaulttextfontsize',fs);
                set(fig,'defaultaxesfontsize',fs);
                set(fig,'defaultuicontrolfontsize',0.9*fs)
                lw = 0.5*fs/10;
                set(fig,'defaultaxeslinewidth',lw)
                set(fig,'defaultlinelinewidth',lw)
                set(fig,'defaultaxesfontname','helvetica')
                set(fig,'defaultaxesfontweight','normal')
                
                set(fig,'user',ud);
                
            case 'export'
                odesdisp = gcf;
                dud = get(odesdisp,'user');
                odset = findobj('name','ODESOLVE Setup');
                sud = get(odset,'user');
                solutions = dud.solutions;
                lsol = length(solutions);
                switch lsol
                    case 0
                        mtitle = 'Solution choice';
                        mmessage = 'There are no solutions.';
                        msgbox(mmessage, mtitle, 'modal')
                        return
                    case 1
                        sol = solutions;
                    otherwise
                        ginput(1);
                        solh = get(odesdisp,'currentobject');
                        if isempty(solh)
                            mtitle = 'Solution choice';
                            mmessage = 'The item selected is not a solution.';
                            msgbox(mmessage, mtitle, 'modal')
                            return
                        end	
                        looking = 1;
                        kk = 0;
                        while looking & kk<=lsol
                            kk = kk + 1;
                            if kk <= lsol
                                kkk = find(solutions(kk).handles == solh);
                                if ~isempty(kkk)
                                    looking = 0;
                                end
                            else
                                looking = 0;
                            end	
                        end
                        if kk <= lsol
                            sol = solutions(kk);
                        else
                            mtitle = 'Solution choice';
                            mmessage = 'The item selected is not a solution.';
                            msgbox(mmessage, mtitle, 'modal')
                            return
                        end
                end
                vars = evalin('base','who');
                no = 1;
                kk = 0;
                while no
                    kk = kk + 1;
                    vstr = ['odedata',num2str(kk)];
                    if ~any(strcmp(vars,vstr))
                        no = 0;
                    end
                end
                tname = sud.c.tname;
                xname = sud.c.xname;
                NN = sud.c.NN;
                ivstr = struct(tname,sol.tdata);
                fieldstr = tname;
                for kk = 1:NN
                    name = xname{kk};
                    name(find(abs(name) == 92)) = [];
                    string = ['ivstr.', name, ' =  sol.xdata(:,kk);'];
                    eval(string)
                    if kk < NN
                        fieldstr = [fieldstr,', ',name];
                    else
                        fieldstr = [fieldstr, ', and ', name];
                    end
                end
                assignin('base',vstr,ivstr);
                mtitle = 'Solution choice';
                mmessage = ['The data has been exported as the sructure ',...
                        vstr, ' with fields ', fieldstr, '.'];
                msgbox(mmessage, mtitle, 'modal')
                
            case 'fileexport'
                
                comp = computer;
                switch  comp
                    case 'PCWIN'
                        prompt = '\*.m';
                    case 'MAC2'
                        prompt = ':*.m';
                    otherwise
                        prompt = ':/.m';
                end
                [fname,pathname] = uiputfile(prompt,'Save the file as:');
                L = length(fname);
                if L > 2 
                    if strcmp(fname([L-1,L]),'.m')
                        fname = fname([1:L-2]);
                    end
                end
                
                [renn,rpnn] = pandecheck('export');
                writefile(fname, pathname, renn, rpnn);
                
            case 'revert'
                
                sud = get(gcf,'user');
                sud.c = sud.o;
                syst = sud.o;
                sh = sud.h;
                set(sh.tvar,'string',syst.tname);
                xname = syst.xname;
                der = syst.der;
                ic = syst.ic;
                NN = syst.NN;
                set(sh.NN,'string',num2str(NN));
                for kk = 1:NN
                    set(sh.xname(kk),'string',xname{kk});
                    set(sh.der(kk),'string',der{kk});
                    set(sh.icval(kk),'string',ic{kk});
                end
                
            case 'quit'
                
                odesolveset = findobj('name','ODESOLVE Setup');
                sud = get(odesolveset,'user');
                if sud.remtd 
                    rmpath(tempdir);
                end
                oldfiles = dir([tempdir,'odestp*.m']);
                for k = 1:length(oldfiles)
                    fn = [tempdir,oldfiles(k).name];
                    fid = fopen(fn,'r');
                    ll = fgetl(fid);
                    ll = fgetl(fid);
                    ll = fgetl(fid);
                    fclose(fid);
                    if strcmp(ll,'%% Created by ODESOLVE')
                        delete(fn)
                    end
                end
                h = findobj('tag','odesolve');
                delete(h);
                
                
            case 'closefcn'
                
                fig = gcf;  
                name = get(fig,'name');
                if strcmp(name,'ODESOLVE Setup') | strcmp(name,'ODESOLVE Display')
                    quest = ['Closing this window will cause all ODESOLVE ',...
                            'windows to close, and ODESOLVE will stop.  ',...
                            'Do you want to quit ODESOLVE?'];
                    butt = questdlg(quest,'Quit ODESOLVE?','Quit','Cancel','Quit');
                    if strcmp(butt,'Quit')
                        odesolve('quit');
                    end
                else
                    delete(findobj('label',name));
                    delete(fig);
                    
                end
                
            case 'outputchoice'
                
                [h,fig] = gcbo;
                value = get(h,'value');
                ud = get(fig,'user');
                NN = ud.c.NN;
                
                
                odeputf = figure('name','ODESOLVE output',...
                    'vis','off',...
                    'numb','off','tag','odesolve');
                odesolve('figdefault',odeputf);
                set(odeputf,'menubar','none');
                switch value
                    case 1     %  Time plot
                        txtstr = ['Check the variables to be plotted vs. ',...
                                ud.c.tname, '.'];
                        pcall = 'odesolve(''tpout'')';
                    case 2    % Phase plane plot
                        txtstr = 'Choose the X-variable and the Y-variable.' ;
                        pcall = 'odesolve(''ppout'')';
                    case 3    % 3 dimensional plot
                        txtstr = 'Choose the X-variable, Y-variable, and  Z-variable.' ;
                        pcall = 'odesolve(''p3out'')';
                end
                canceldata = ud.h.output;
                
                figdir = uicontrol('style','text',...
                    'horiz','left',...
                    'string',txtstr);
                ext = get(figdir,'extent');
                nudge = 2;
                left = 2;
                bbott = 2;
                sep = 2;
                ht = ext(4)+nudge;
                fdw = ext(3) + nudge;
                figw = fdw +2*left;
                buttw = fdw/2;
                cbutt = uicontrol('style','push',...
                    'pos',[left,bbott,buttw,ht],...
                    'string','Cancel',...
                    'call','odesolve(''outcancel'')',...
                    'user',canceldata);
                pbutt = uicontrol('style','push',...
                    'pos',[left+buttw,bbott,buttw,ht],...
                    'string','Change',...
                    'user',fig,...
                    'call',pcall);
                switch value
                    case 1     %  Time series
                        checkbott = bbott + ht + sep;
                        checkleft = left + 10;
                        checkw = fdw - checkleft;
                        cpos = [checkleft,checkbott,checkw,ht];
                        allcall = [
                            '[h,fig] = gcbo;',...
                                'checks = get(fig,''user'');',...
                                'NN = length(checks)-1;',...
                                'allch = checks(NN+1);',...
                                'val = get(allch,''value'');',...
                                'set(checks(1:NN),''value'',val);'];
                        checks(NN+1) = uicontrol('style','check',...
                            'pos',cpos,...
                            'string','all',...
                            'value',1,...
                            'call',allcall);
                        checkbott = checkbott +ht+sep;
                        ichcall = [
                            '[h,fig] = gcbo;',...
                                'checks = get(fig,''user'');',...
                                'NN = length(checks)-1;',...
                                'allch = checks(NN+1);',...
                                'values = get(checks(1:NN),''value'');',...
                                'vals = zeros(1:NN);',...
                                'for kk = 1:NN,',...
                                '  vals(kk) = values{kk};',...
                                'end,',...
                                'if all(vals),',...
                                '  set(allch,''value'',1),',...
                                'else,',...
                                '  set(allch,''value'',0),',...
                                'end'];
                        
                        
                        for k = 1:NN
                            cpos = [checkleft,checkbott,checkw,ht];
                            checks(NN+1-k) = uicontrol('style','check',...
                                'pos',cpos,...
                                'value',1,...
                                'string',ud.c.xname{NN+1-k},...
                                'call',ichcall);
                            checkbott = checkbott+ht;
                        end
                        fdbott = checkbott+sep;
                        putdata = checks;
                        
                    case 2  % Phase plane plot
                        radbott = bbott + ht + sep;
                        xradleft = left +10;
                        radw = fdw/2 - 10;
                        yradleft = xradleft +radw+10;
                        for k = 1:NN
                            xradpos = [xradleft,radbott,radw,ht];
                            yradpos = [yradleft,radbott,radw,ht];
                            if k == NN;
                                xval = 1;
                                yval = 0;
                            elseif k == NN - 1;
                                xval = 0;
                                yval = 2;
                            else
                                xval = 0;
                                yval = 0;
                            end
                            xrad(NN+1-k) = uicontrol('style','radio',...
                                'pos',xradpos,...
                                'string',ud.c.xname{NN+1-k},...
                                'value',xval,...
                                'max',NN+1-k);
                            yrad(NN+1-k) = uicontrol('style','radio',...
                                'pos',yradpos,...
                                'string',ud.c.xname{NN+1-k},...
                                'value',yval,...
                                'max',NN+1-k);
                            radbott = radbott + ht;
                        end
                        for k = 1:NN
                            set(xrad(k),'user',[xrad(:,[1:(k-1),(k+1):NN]),yrad(k)]);
                            set(yrad(k),'user',[yrad(:,[1:(k-1),(k+1):NN]),xrad(k)]);
                        end
                        rcall = [
                            'me = get(gcf,''currentobject'');',...
                                'kk = get(me,''max'');',...
                                'set(get(me,''user''),''value'',0),',...
                                'set(me,''value'',kk);'];
                        set([xrad,yrad],'call',rcall);
                        xpos = [xradleft,radbott,radw,ht];
                        ypos = [yradleft,radbott,radw,ht];
                        uicontrol('style','text',...
                            'pos',xpos,...
                            'string','X-variable');
                        uicontrol('style','text',...
                            'pos',ypos,...
                            'string','Y-variable');
                        radbott = radbott +ht;
                        fdbott = radbott+sep;
                        putdata = struct('xvar',xrad,'yvar',yrad);
                        
                    case 3  % 3 dimensional plot
                        radbott = bbott + ht + sep;
                        xradleft = left +10;
                        radw = fdw/3 - 10;
                        yradleft = xradleft +radw+10;
                        zradleft = yradleft +radw+10;
                        if NN == 1
                            errordlg(['A 3 dimensional plot is not useful for a single' ...
                                    ' equation.'])
                            return;
                        end
                        xradpos = [xradleft,radbott,radw,ht];
                        yradpos = [yradleft,radbott,radw,ht];
                        zradpos = [zradleft,radbott,radw,ht];
                        radbott = radbott +ht;
                        xrad(NN+1) = uicontrol('style','radio',...
                            'pos',xradpos,...
                            'string',ud.c.tname,...
                            'value',0,...
                            'max',NN+1);
                        yrad(NN+1) = uicontrol('style','radio',...
                            'pos',yradpos,...
                            'string',ud.c.tname,...
                            'value',0,...
                            'max',NN+1);
                        zrad(NN+1) = uicontrol('style','radio',...
                            'pos',zradpos,...
                            'string',ud.c.tname,...
                            'value',NN+1,...
                            'max',NN+1);
                        
                        for k = 1:NN
                            xradpos = [xradleft,radbott,radw,ht];
                            yradpos = [yradleft,radbott,radw,ht];
                            zradpos = [zradleft,radbott,radw,ht];
                            if k == NN;
                                xval = 1;
                                yval = 0;
                                zval = 0;
                            elseif k == NN - 1;
                                xval = 0;
                                yval = 2;
                                zval = 0;
                            else
                                xval = 0;
                                yval = 0;
                                zval = 0;
                            end
                            xrad(NN+1-k) = uicontrol('style','radio',...
                                'pos',xradpos,...
                                'string',ud.c.xname{NN+1-k},...
                                'value',xval,...
                                'max',NN+1-k);
                            yrad(NN+1-k) = uicontrol('style','radio',...
                                'pos',yradpos,...
                                'string',ud.c.xname{NN+1-k},...
                                'value',yval,...
                                'max',NN+1-k);
                            zrad(NN+1-k) = uicontrol('style','radio',...
                                'pos',zradpos,...
                                'string',ud.c.xname{NN+1-k},...
                                'value',zval,...
                                'max',NN+1-k);
                            radbott = radbott + ht;
                        end
                        
                        for k = 1:NN+1
                            set(xrad(k),'user',[xrad(:,[1:(k-1),(k+1):NN+1]),yrad(k),zrad(k)]);
                            set(yrad(k),'user',[yrad(:,[1:(k-1),(k+1):NN+1]),xrad(k),zrad(k)]);
                            set(zrad(k),'user',[zrad(:,[1:(k-1),(k+1):NN+1]),xrad(k),yrad(k)]);
                        end
                        rcall = [
                            'me = get(gcf,''currentobject'');',...
                                'kk = get(me,''max'');',...
                                'set(get(me,''user''),''value'',0),',...
                                'set(me,''value'',kk);'];
                        set([xrad,yrad,zrad],'call',rcall);
                        xpos = [xradleft-8,radbott,radw,ht];
                        ypos = [yradleft-8,radbott,radw,ht];
                        zpos = [zradleft-8,radbott,radw,ht];
                        uicontrol('style','text',...
                            'pos',xpos,...
                            'string','X-variable');
                        uicontrol('style','text',...
                            'pos',ypos,...
                            'string','Y-variable');
                        uicontrol('style','text',...
                            'pos',zpos,...
                            'string','Z-variable');
                        radbott = radbott +ht;
                        fdbott = radbott+sep;
                        putdata = struct('xvar',xrad,'yvar',yrad,'zvar',zrad);
                        
                end  % switch value
                
                set(figdir,'pos',[left,fdbott,fdw,ht]);
                fight = fdbott + ht + nudge;
                set(odeputf,'pos',[100,300,figw,fight],...
                    'user',putdata,...
                    'vis','on');
                set(get(odeputf,'child'),'vis','on');
                set(findobj(odeputf,'type','uicontrol'),'units','normal');
                figure(odeputf)
                
                
            case 'outcancel'  % Callback for Cancel in output choice window.
                
                [h,fig] = gcbo;
                cud =get(h,'user');
                set(cud,'value',get(cud,'user'));
                close;
                
                
                
                
            case 'tpout'  % setup time plot output
                
                [h,fig] = gcbo;
                checks = get(fig,'user');
                NN = length(checks) - 1;
                checks = checks(1:NN);
                val = get(checks,'value');
                val = c2m(val);
                NNN = find(val);
                if isempty(NNN)
                    errordlg('At least one variable must be chosen.');
                    return;
                end
                odesetup = get(h,'user');
                ud = get(odesetup,'user');
                ud.NNN = NNN;     
                hh = ud.h.outputstr;
                if length(NNN) == ud.c.NN
                    str = ['     all variables vs ', ud.c.tname,'.'];
                elseif NNN == 1
                    str = ['     ', ud.c.xname{NNN}, ' vs ' ud.c.tname,'.'];
                else
                    str = ['     ', ud.c.xname{NNN(1)}];
                    for kk = 2:length(NNN)
                        str = [str, ', ', ud.c.xname{NNN(kk)}];
                    end
                    str = [str, ' vs ' ud.c.tname,'.'];
                end
                set(hh,'string',str);     
                ud.outputhandle = @tpout;
                set(ud.h.sobutt,'enable','off');
                set(ud.h.somenu,'enable','off');
                set(odesetup,'user',ud);
                close
                
            case 'ppout'
                
                [h,fig] = gcbo;
                ppud = get(fig,'user');
                xval = c2m(get(ppud.xvar,'value'));
                yval = c2m(get(ppud.yvar,'value'));
                xvar = find(xval);
                yvar = find(yval);
                if isempty(xvar) | isempty(yvar)
                    errordlg('You must choose both an X-variable and a Y-variable.');
                    return;
                end
                odesetup = get(h,'user');
                ud = get(odesetup,'user');
                hh = ud.h.outputstr;
                str = ['     ', ud.c.xname{yvar}, ' vs ',ud.c.xname{xvar}];
                set(hh,'string',str);
                ud.outputhandle = @ppout;
                set(ud.h.sobutt,'enable','off');
                set(ud.h.somenu,'enable','off');
                %  ud.outputhandle = @odesp2;
                ud.NNN = [xvar,yvar];
                set(odesetup,'user',ud);
                close
                
            case 'p3out'
                
                [h,fig] = gcbo;
                ppud = get(fig,'user');
                xval = c2m(get(ppud.xvar,'value'));
                yval = c2m(get(ppud.yvar,'value'));
                zval = c2m(get(ppud.zvar,'value'));
                xvar = find(xval);
                yvar = find(yval);
                zvar = find(zval);
                if isempty(xvar) | isempty(yvar) | isempty(zvar)
                    errordlg('You must choose an X-variable,a Y-variable, and a Z-variable.');
                    return;
                end
                odesetup = get(h,'user');
                ud = get(odesetup,'user');
                hh = ud.h.outputstr;
                NN = ud.c.NN;
                if xvar == NN+1
                    xn = ud.c.tname;
                else 
                    xn = ud.c.xname{xvar};
                end
                if yvar == NN+1
                    yn = ud.c.tname;
                else 
                    yn = ud.c.xname{yvar};
                end
                if zvar == NN+1
                    zn = ud.c.tname;
                else 
                    zn = ud.c.xname{zvar};
                end
                
                str = ['     ', xn, ', ', yn, ', and ', zn];
                set(hh,'string',str);
                ud.outputhandle = @p3out;
                set(ud.h.sobutt,'enable','off');
                set(ud.h.somenu,'enable','off');
                ud.NNN = [xvar,yvar,zvar];
                set(odesetup,'user',ud);
                close
                
            case 'settings'
                
                % Routine for changing solver settings.
                
                [h,fig] = gcbo;
                sud = get(fig,'user');
                NN = sud.c.NN;
                settings = sud.settings;
                ecol = sud.ecolor;
                setf = findobj('name','Solver settings');
                if ~isempty(setf)
                    close(setf)
                end
                setf = figure('name','Solver settings',...
                    'vis','off',...
                    'tag','odesolve',...
                    'numb','off');
                odesolve('figdefault',setf);
                set(setf,'menubar','none')
                val = get(sud.h.solver,'value');
                solver = sud.solvers{val};
                titlestring = ['Settings for ', solver,'.'];
                title = uicontrol('style','text',...
                    'horizon','center',...
                    'string',['Settings for ', solver, '.'],...
                    'visible','off');
                ext = get(title,'extent');
                delete(title)
                texth = ext(4) + 4;
                left = 1;
                eleft = left + 4;
                frw = 258;
                figw = frw + 2*left;
                frsep = 1;
                sep = 1;
                
                %% The buttons.
                
                buttw = (frw-3*frsep)/2;
                cwind = [left+frsep,1,buttw,texth];
                awind = [cwind(1)+cwind(3)+frsep,1,buttw,texth];
                
                ccall = [
                    '[h,fig] = gcbo;',...
                        'close(fig)'];
                butt(1) = uicontrol('style','push',...
                    'pos',cwind,...
                    'string','Cancel',...
                    'call',ccall,...
                    'visible','off');
                butt(2) = uicontrol('style','push',...
                    'pos',awind,...
                    'string','Apply',...
                    'call','odesolve(''setapply'')',...
                    'visible','off');
                
                frbot = 1 + texth + frsep;
                KK = ceil(NN/2);
                if strcmp(solver,'ode15s')
                    frht = (9 + KK)*(texth+sep)+3*sep;
                else
                    frht = (7 + KK)*(texth+sep)+sep;
                end
                frwind = [left,frbot,frw,frht];
                fr = uicontrol('style','frame',...
                    'pos', frwind);
                fight = frbot + frht + frsep;
                figwind = [20,200,figw,fight];
                set(setf,'pos',figwind);
                
                optw = 170;
                ncbot = frbot+sep;
                ncwind = [eleft,ncbot,optw,texth];
                normcont = uicontrol('style','text',...
                    'horizon','left',...
                    'string','Norm control',...
                    'visible','off',...
                    'pos',ncwind);
                editleft = eleft + optw;
                editw = 80;
                encwind = [editleft,ncbot+2,editw,texth];
                onoffstr ={'       on','      off'};
                editnormcont = uicontrol('style','popupmenu',...
                    'horizon','center',...
                    'string',onoffstr,...
                    'value',settings.ncval,...
                    'visible','off',...
                    'backgroundcolor',ecol,...
                    'pos',encwind);
                
                mscall = [
                    '[h,fig] = gcbo;',...
                        'str = get(h,''string'');',...
                        'nstr = str2num(str);',...
                        'if isempty(nstr) | nstr <= 0,',...
                        '  set(h,''string'',''computed'');'...
                        'end,',...
                    ];
                
                msbot = ncbot + texth + sep;
                mswind = [eleft,msbot,optw,texth];
                maxstep = uicontrol('style','text',...
                    'horizon','left',...
                    'string','Maximum step size',...
                    'visible','off',...
                    'pos',mswind);
                emswind = [editleft,msbot+2,editw,texth];
                editmaxstep = uicontrol('style','edit',...
                    'horizon','center',...
                    'string',settings.mstep,...
                    'call',mscall,...
                    'visible','off',...
                    'backgroundcolor',ecol,...
                    'pos',emswind);
                iscall = [
                    '[h,fig] = gcbo;',...
                        'str = get(h,''string'');',...
                        'nstr = str2num(str);',...
                        'if isempty(nstr) | nstr <= 0,',...
                        '  set(h,''string'',''computed'');'...
                        'end,',...
                    ];
                isbot = msbot+texth+sep;
                iswind = [eleft,isbot,optw,texth];
                initstep = uicontrol('style','text',...
                    'horizon','left',...
                    'string','Initial step size',...
                    'visible','off',...
                    'pos',iswind);
                eiswind = [editleft,isbot+2,editw,texth];
                editinitstep = uicontrol('style','edit',...
                    'horizon','center',...
                    'string',settings.istep,...
                    'call',iscall,...
                    'visible','off',...
                    'backgroundcolor',ecol,...
                    'pos',eiswind);
                refcall = [
                    '[h,fig] = gcbo;',...
                        'refine = str2num(get(h,''string''));',...
                        'if isempty(refine) | refine<=0,',...
                        '  refine = 1;',... 
                        'end,',...
                        'refine = ceil(refine);',...
                        'set(h,''string'',refine);',...
                    ];
                
                refbot = isbot + texth + sep;
                refwind = [eleft,refbot,optw,texth];
                refine = uicontrol('style','text',...
                    'horizon','left',...
                    'string','Refine',...
                    'visible','off',...
                    'pos',refwind);
                erefwind = [editleft,refbot+2,editw,texth];
                editrefine = uicontrol('style','edit',...
                    'horizon','center',...
                    'string',settings.refine,...
                    'call',refcall,...
                    'visible','off',...
                    'backgroundcolor',ecol,...
                    'pos',erefwind);
                rtolbot = refbot + texth + sep;
                rtolwind = [eleft,rtolbot,optw, texth];
                rtol = uicontrol('style','text',...
                    'horizon','left',...
                    'string','Relative tolerance',...
                    'visible','off',...
                    'pos',rtolwind);
                ertolwind = [editleft,rtolbot+2,editw,texth];
                editrtol = uicontrol('style','edit',...
                    'horizon','center',...
                    'string',num2str(settings.rtol),...
                    'visible','off',...
                    'backgroundcolor',ecol,...
                    'pos',ertolwind);
                
                vatolcall = [
                    '[h,fig] = gcbo;',...
                        'ud = get(fig,''user'');',...
                        'set(ud.editatol,''string'','''');',...
                    ];
                
                vleft = eleft +8;
                vw = 40;
                evw = 80;
                vjbot = rtolbot;
                xname = sud.c.xname;
                for kk = 1:KK
                    vjbot = vjbot +texth + sep;
                    for jj = 0:1
                        vk = 2*(KK - kk) + 1 + jj;
                        if vk <= NN
                            vjleft = vleft + jj*(vw + evw + sep);
                            vjwind = [vjleft,vjbot,vw,texth];
                            tv = uicontrol('style','text',...
                                'horizon','center',...
                                'string',xname{vk},...
                                'visible','off',...
                                'pos',vjwind);
                            evjleft = vjleft + vw + sep;
                            evjwind = [evjleft,vjbot+2,evw,texth];
                            ev(vk) = uicontrol('style','edit',...
                                'horizon','center',...
                                'string','???',...
                                'call',vatolcall,...
                                'visible','off',...
                                'backgroundcolor',ecol,...
                                'pos',evjwind);
                        end
                    end
                end
                atolbot = vjbot + texth + sep;
                atolwind = [eleft,atolbot,optw, texth];
                atol = uicontrol('style','text',...
                    'horizon','left',...
                    'string','Absolute tolerance:              All',...
                    'visible','off',...
                    'pos',atolwind);
                eatolwind = [editleft,atolbot+2,editw,texth];
                
                atolcall = [
                    '[h,fig] = gcbo;',...
                        'ud = get(fig,''user'');',...
                        'ev = ud.ev;',...
                        'at = get(ud.editatol,''string'');',...	'[h,fig] = gcbo;',...
                        'ud = get(fig,''user'');',...
                        'ev = ud.ev;',...
                        'at = get(ud.editatol,''string'');',...
                        'set(ev,''string'',at);',...
                    ];
                
                
                editatol = uicontrol('style','edit',...
                    'horizon','center',...
                    'string','what goes here?',...
                    'visible','off',...
                    'call',atolcall,...
                    'backgroundcolor',ecol,...
                    'pos',eatolwind);
                atol = settings.atol;
                if length(atol) == 1
                    set(ev,'string',num2str(atol));
                    set(editatol,'string',atol);
                else
                    for kk = 1:NN
                        set(ev(kk),'string',num2str(atol(kk)));
                        set(editatol,'string','');
                    end
                end
                
                tbot = atolbot + texth + sep;
                tw = figw - 2*eleft;
                twind = [eleft,tbot,tw,texth];
                title = uicontrol('style','text',...
                    'horizon','center',...
                    'string',['Settings for ', solver, '.'],...
                    'visible','off',...
                    'pos',twind);
                
                ud.editnormcont = editnormcont;
                ud.editmaxstep = editmaxstep;
                ud.editinitstep = editinitstep;
                ud.editrefine = editrefine;
                ud.editrtol = editrtol;
                ud.ev = ev;
                ud.editatol = editatol;
                set(setf,'user',ud);
                
                
                
                
                set(setf,'vis','on');
                set(get(setf,'child'),'vis','on');
                
            case 'setapply'
                
                [h,fig] = gcbo;
                ud = get(fig,'user');
                odset = findobj('name','ODESOLVE Setup');
                sud = get(odset,'user');
                settings = sud.settings;
                atol = get(ud.editatol,'string');
                if isempty(atol);
                    ev = ud.ev;
                    atol = zeros(1,sud.c.NN);
                    for kk = 1:sud.c.NN
                        atol(kk) = str2num(get(ev(kk),'string'));
                    end
                else
                    atol = str2num(atol);
                end
                settings.atol = atol;
                str = get(ud.editrtol,'string');
                settings.rtol = str2num(get(ud.editrtol,'string'));
                settings.refine = str2num(get(ud.editrefine,'string'));
                is = get(ud.editinitstep,'string');
                if strcmp(is,'computed')
                    settings.istep = is;
                else
                    settings.istep = str2num(is);
                end
                ms = get(ud.editmaxstep,'string');
                if strcmp(ms,'computed')
                    settings.mstep = ms;
                else
                    settings.mstep = str2num(ms);
                end
                settings.ncval = get(ud.editnormcont,'value');
                sud.settings = settings;
                set(odset,'user',sud);
                sud = get(odset,'user');
                close(fig);
                
            case 'restart'
                
                odesolve
                
            end
            
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            function status = tpout(t,y,flag,varargin)
            
            status = 0;                             % Assume stop button wasn't pushed.
            chunk = 128;                            % Memory is allocated in chunks.
            odset = findobj('name','ODESOLVE Setup');
            sud = get(odset,'user');
            NNN = sud.NNN;     % The variables to be plotted.
            dispfig = sud.h.disp;
            ud = get(dispfig,'user');
            NS = length(ud.solutions);  % The number of solutions already plotted.
            
            if nargin < 3 | isempty(flag) % odesplot(t,y) [v5 syntax] or odesplot(t,y,'')
                
                % Append t and y to ud.t and ud.y, allocating new memory if necessary.
                nt = length(t);
                chunk = max(chunk,nt);
                [rows,cols] = size(ud.y);
                plsols = length(NNN);   % The number of plots 
                oldi = ud.i;    % The number of previous steps
                newi = oldi + nt;    % The current number of steps
                if newi > rows
                    ud.t = [ud.t; zeros(chunk,1)];
                    ud.y = [ud.y; zeros(chunk,cols)];
                end
                ud.t(oldi+1:newi) = t;  
                ud.y(oldi+1:newi,:) = y.';
                ud.i = newi;  % The current number of steps
                set(dispfig,'UserData',ud);
                
                if ud.stop == 1                       % Has stop button been pushed?
                    status = 1;
                    tpout([],[],'done');
                else
                    % Rather than redraw all of the data every timestep, we will simply move
                    % the line segments for the new data, not erasing.  But if the data has
                    % moved out of the axis range, we redraw everything.  If there is
                    % already a solution plotted, we only move the line seqments.
                    
                    ylim = get(gca,'ylim');
                    if (NS == 0) & ((oldi == 1) | (min(y(:)) < ylim(1)) | (ylim(2) < max(y(:))))
                        % Replot everything if this is the first plot, and the data is
                        % out of axis range or if the plot is just initialized. 
                        for j = 1:plsols
                            set(ud.lines(j),'Xdata',ud.t(1:newi),'Ydata',ud.y(1:newi,NNN(j)));
                        end
                    else
                        % Plot only the new data if the axis range is ok, or if this is
                        % plotted on top of another solution.
                        for j = 1:plsols
                            set(ud.line(j),'Xdata',ud.t(oldi:newi),'Ydata',ud.y(oldi:newi,NNN(j)));
                        end
                    end
                end
                
            else
                switch(flag)
                    case 'init'                           % odesplot(tspan,y0,'init')
                        cols = length(y);
                        ud.t = zeros(chunk,1);
                        ud.y = zeros(chunk,cols);
                        ud.i = 1;
                        ud.t(1) = t(1);
                        ud.y(1,:) = y.';
                        if sud.pplotflag
                            mm = '.';
                        else
                            mm = 'none';
                        end
                        
                        axes(ud.axes);
                        
                        % Rather than redraw all data at every timestep, we will simply move
                        % the last line segment along, not erasing it.
                        
                        ms = 6;  % Marker size
                        ud.lines = plot(ud.t(1),ud.y(1,NNN),'-',...
                            'marker',mm,...
                            'markersize',ms,...
                            'erasemode','normal');
                        hold on
                        
                        ud.line = plot(ud.t(1),ud.y(1,NNN),'-',...
                            'marker',mm,...
                            'markersize',ms,...
                            'EraseMode','none');
                        set(gca,'XLim',sud.c.tvals(2:3));
                        hold off
                        
                        % The STOP button.
                        h = findobj(dispfig,'Tag','stop');
                        if isempty(h)
                            ud.stop = 0;
                            pos = get(0,'DefaultUicontrolPosition');
                            pos(1) = pos(1) - 15;
                            pos(2) = pos(2) - 15;
                            str = [
                                '[h,dispfig] = gcbo;',...
                                    'ud=get(dispfig,''UserData''); ',...
                                    'ud.stop=1;',...
                                    'set(dispfig,''UserData'',ud);'];
                            uicontrol( ...
                                'Style','push', ...
                                'String','Stop', ...
                                'Position',pos, ...
                                'Callback',str, ...
                                'Tag','stop');
                        else
                            set(h,'Visible','on');            % make sure it's visible
                            ud.stop = 0;
                        end
                        set(dispfig,'UserData',ud);
                        
                    case 'done'                           
                        ud.t = ud.t(1:ud.i);
                        ud.y = ud.y(1:ud.i,:);
                        delete(ud.line);
                        ud.line=[];
                        for j = 1:length(NNN)
                            set(ud.lines(j),'xdata',ud.t,...
                                'ydata',ud.y(:,NNN(j)),...
                                'erase','normal');
                        end
                        set(dispfig,'UserData',ud);
                        set(findobj(dispfig,'Tag','stop'),'Visible','off');
                        % refresh;                          % redraw figure to remove marker frags
                end
            end
            
            drawnow;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            function status = ppout(t,y,flag,varargin)
            
            status = 0;                             % Assume stop button wasn't pushed.
            chunk = 128;                            % Memory is allocated in chunks.
            odset = findobj('name','ODESOLVE Setup');
            sud = get(odset,'user');
            kk1 = sud.NNN(1);     % The variables to be plotted.
            kk2 = sud.NNN(2);     % The variables to be plotted.
            dispfig = sud.h.disp;
            ud = get(dispfig,'user');
            NS = length(ud.solutions);  % The number of solutions already plotted.
            
            if nargin < 3 | isempty(flag) % odesp2(t,y) [v5 syntax] or odesp2(t,y,'')
                
                % Append y to ud.y, allocating new memory if necessary.
                nt = length(t);
                chunk = max(chunk,nt);
                rows = size(ud.y,1);
                oldi = ud.i;
                newi = oldi + nt;
                if newi > rows
                    ud.y = [ud.y; zeros(chunk,2)];
                end
                y = y([kk1 kk2],:);  % Keep only the two rows of interest.
                ud.y(oldi+1:newi,:) = y.';
                ud.i = newi;
                set(dispfig,'UserData',ud);
                
                if ud.stop == 1                       % Has stop button been pushed?
                    status = 1;
                    ppout([],[],'done');
                else
                    % Rather than redraw all of the data every timestep, we will simply move
                    % the line segments for the new data, not erasing.  But if the data has
                    % moved out of the axis range, we redraw everything.  If there is
                    % already a solution plotted, we only move the line seqments.
                    
                    xlim = get(gca,'xlim');
                    ylim = get(gca,'ylim');
                    % Replot everything if out of axis range or if just initialized.
                    if ~NS & ((oldi == 1) | ...
                            (min(y(1,:)) < xlim(1)) | (xlim(2) < max(y(1,:))) | ...
                            (min(y(2,:)) < ylim(1)) | (ylim(2) < max(y(2,:))))
                        % Replot everything if this is the first plot, and the data is
                        % out of axis range or if the plot is just initialized. 
                        set(ud.lines, ...
                            'Xdata',ud.y(1:newi,1), ...
                            'Ydata',ud.y(1:newi,2));
                    else
                        % Plot only the new data if the axis range is ok, or if this is
                        % plotted on top of another solution.
                        
                        % set(ud.line,'Color',co(1,:));     % "erase" old segment
                        set(ud.line, ...
                            'Xdata',ud.y(oldi:newi,1), ...
                            'Ydata',ud.y(oldi:newi,2));
                    end
                end
                
            else
                switch(flag)
                    case 'init'                           % odesp2(tspan,y0,'init')
                        ud.y = zeros(chunk,2);
                        ud.i = 1;
                        ud.y(1,:) = y([kk1 kk2]).';
                        if sud.pplotflag
                            mm = '.';
                        else
                            mm = 'none';
                        end
                        axes(ud.axes);
                        
                        % Rather than redraw all data at every timestep, we will simply move
                        % the last line segment along, not erasing it.
                        ms = 6;  % Marker size
                        
                        ud.lines = plot(y(kk1),y(kk2),'-',...
                            'marker',mm,...
                            'markersize',ms);
                        hold on
                        ud.line = plot(y(kk1),y(kk2),'-',...
                            'marker',mm,...
                            'markersize',ms,...
                            'EraseMode','none');
                        hold off
                        
                        % The STOP button.
                        h = findobj(dispfig,'Tag','stop');
                        if isempty(h)
                            ud.stop = 0;
                            pos = get(0,'DefaultUicontrolPosition');
                            pos(1) = pos(1) - 15;
                            pos(2) = pos(2) - 15;
                            str = [
                                '[h,dispfig] = gcbo;',...
                                    'ud=get(dispfig,''UserData''); ',...
                                    'ud.stop=1;',...
                                    'set(dispfig,''UserData'',ud);'];
                            
                            uicontrol( ...
                                'Style','push', ...
                                'String','Stop', ...
                                'Position',pos, ...
                                'Callback',str, ...
                                'Tag','stop');
                        else
                            set(h,'Visible','on');            % make sure it's visible
                            ud.stop = 0;
                        end
                        set(dispfig,'UserData',ud);
                        
                    case 'done'                           % odesp2([],[],'done')
                        ud.y = ud.y(1:ud.i,:);
                        set(ud.lines, ...
                            'Xdata',ud.y(:,1), ...
                            'Ydata',ud.y(:,2));
                        set(dispfig,'UserData',ud);
                        set(findobj(dispfig,'Tag','stop'),'Visible','off');
                        refresh;                          % redraw figure to remove marker frags
                        
                end
            end
            
            drawnow;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            function status = p3out(t,y,flag,varargin)
            
            status = 0;                             % Assume stop button wasn't pushed.
            chunk = 128;                            % Memory is allocated in chunks.
            odset = findobj('name','ODESOLVE Setup');
            sud = get(odset,'user');
            NNN = sud.NNN;     % The variables to be plotted.
            NN = sud.c.NN;    % The number of x-variables.
            dispfig = sud.h.disp;
            ud = get(dispfig,'user');
            NS = length(ud.solutions);  % The number of solutions already plotted.
            
            if nargin < 3 | isempty(flag) % odesp2(t,y) [v5 syntax] or odesp2(t,y,'')
                
                % Append appropriate parts of y and t to ud.y, allocating new memory
                % if necessary. 
                nt = length(t);
                chunk = max(chunk,nt);
                rows = size(ud.y,1);
                oldi = ud.i;
                newi = oldi + nt;
                for k = 1:3  % Keep only the rows of interest.
                    if NNN(k) > NN;
                        yy(:,k) = t;
                    else
                        yy(:,k) = y(NNN(k),:).';
                    end
                end
                if newi > rows
                    ud.y = [ud.y; zeros(chunk,3)];
                end
                
                ud.y(oldi+1:newi,:) = yy;
                ud.i = newi;
                set(dispfig,'UserData',ud);
                
                if ud.stop == 1                       % Has stop button been pushed?
                    status = 1;
                    p3out([],[],'done');
                else
                    % Rather than redraw all of the data every timestep, we will simply move
                    % the line segments for the new data, not erasing.  But if the data has
                    % moved out of the axis range, we redraw everything.  If there is
                    % already a solution plotted, we only move the line seqments.
                    
                    xlim = get(gca,'xlim');
                    ylim = get(gca,'ylim');
                    zlim = get(gca,'zlim');
                    % Replot everything if out of axis range or if just initialized.
                    if ~NS & ((oldi == 1) | ...
                            (min(yy(:,1)) < xlim(1)) | (xlim(2) < max(yy(:,1))) | ...
                            (min(yy(:,2)) < ylim(1)) | (ylim(2) < max(yy(:,2))) | ...
                            (min(yy(:,3)) < zlim(1)) | (zlim(2) < max(yy(:,3))))
                        % Replot everything if this is the first plot, and the data is
                        % out of axis range or if the plot is just initialized. 
                        set(ud.lines, ...
                            'Xdata',ud.y(1:newi,1), ...
                            'Ydata',ud.y(1:newi,2),...
                            'Zdata',ud.y(1:newi,3));
                    else
                        % Plot only the new data if the axis range is ok, or if this is
                        % plotted on top of another solution.
                        
                        % set(ud.line,'Color',co(1,:));     % "erase" old segment
                        set(ud.line, ...
                            'Xdata',ud.y(oldi:newi,1), ...
                            'Ydata',ud.y(oldi:newi,2), ...
                            'Zdata',ud.y(oldi:newi,3));
                    end
                end
                
            else
                switch(flag)
                    
                    case 'init'                           % odesp2(tspan,y0,'init')
                        ud.y = zeros(chunk,3);
                        ud.i = 1;
                        for k = 1:3  % Keep only the rows of interest.
                            if NNN(k) > NN;
                                yy(:,k) = t(1);
                            else
                                yy(:,k) = y(NNN(k),:).';
                            end
                        end
                        ud.y(1,:) = yy;
                        if sud.pplotflag
                            mm = '.';
                        else
                            mm = 'none';
                        end
                        axes(ud.axes);
                        
                        % Rather than redraw all data at every timestep, we will simply move
                        % the last line segment along, not erasing it.
                        ms = 6;  % Marker size
                        
                        ud.lines = plot3(yy(1),yy(2),yy(3),'-',...
                            'marker',mm,...
                            'markersize',ms);
                        hold on
                        ud.line = plot3(yy(1),yy(2),yy(3),'-',...
                            'marker',mm,...
                            'markersize',ms,...
                            'EraseMode','none');
                        hold off
                        
                        % The STOP button.
                        h = findobj(dispfig,'Tag','stop');
                        if isempty(h)
                            ud.stop = 0;
                            pos = get(0,'DefaultUicontrolPosition');
                            pos(1) = pos(1) - 15;
                            pos(2) = pos(2) - 15;
                            str = [
                                '[h,dispfig] = gcbo;',...
                                    'ud=get(dispfig,''UserData''); ',...
                                    'ud.stop=1;',...
                                    'set(dispfig,''UserData'',ud);'];
                            
                            uicontrol( ...
                                'Style','push', ...
                                'String','Stop', ...
                                'Position',pos, ...
                                'Callback',str, ...
                                'Tag','stop');
                        else
                            set(h,'Visible','on');            % make sure it's visible
                            ud.stop = 0;
                        end
                        set(dispfig,'UserData',ud);
                        
                    case 'done'                           % odesp2([],[],'done')
                        ud.y = ud.y(1:ud.i,:);
                        set(ud.lines, ...
                            'Xdata',ud.y(:,1), ...
                            'Ydata',ud.y(:,2), ...
                            'Zdata',ud.y(:,3));
                        set(dispfig,'UserData',ud);
                        set(findobj(dispfig,'Tag','stop'),'Visible','off');
                        refresh;                          % redraw figure to remove marker frags
                        
                end
            end
            
            drawnow;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            function M = c2m(c)
            
            [rows,cols] = size(c);
            if rows*cols == 0
                M = [];
            else
                M = zeros(rows,cols);
                for i=1:rows 
                    M(i,:) = [c{i,[1:cols]}]; 
                end
            end
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            function outstring = pareval(string,par,value)
            
            string = deblank(string);
            string(find(abs(string)==32))=[];
            par = deblank(par);
            value = ['(',value,')'];
            ll = length(string);
            lp = length(par);
            lv = length(value);
            if strcmp(par,string)
                string = value;
            elseif (ll >= lp+1)
                k = findstr(par,string);
                lk = length(k);
                lopstr = '(+-*/^';
                ropstr = ')+-*/^';
                s = [];
                pos = 1;
                for jk = 1:lk
                    if (((k(jk) == 1)|(find(lopstr == string(k(jk)-1))))...
                            &((k(jk)+lp-1 == ll)|(find(ropstr == string(k(jk) + lp)))))
                        s = [s,string(pos:(k(jk)-1)),value];
                        pos = k(jk)+lp;
                    end
                end
                string = [s,string(pos:ll)];
            end
            outstring = string;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            function [renn,rpnn] = pandecheck(flag)
            
            %% Checks which parameters and expressions are actually used.
            
            odesolveset = findobj('name','ODESOLVE Setup');
            sud = get(odesolveset,'user');
            
            %% The expressions.  Expressions can contain parameters.  Find the
            %% expressions that are actually used in a derivative.
            
            exname = sud.h.exname;
            exval = sud.h.exval;
            der = sud.c.der;
            NN = sud.c.NN;
            renn = [];
            for kk = 1:length(exname)
                exn = char(get(exname(kk),'string'));
                if ~isempty(exn)
                    present = 0;
                    kkk = 1;
                    while present == 0 & kkk <= NN
                        if ~isempty(findstr(exn,char(der(kkk))))
                            present = 1;
                        end
                        kkk = kkk + 1;
                    end
                    if present
                        renn = [renn,kk];
                    end
                end
            end
            krex = length(renn);
            
            %% Set up the parameters.  First find those that are actually used in
            %% either an expression or a drivative.
            
            pname = sud.h.pname;
            icond = sud.c.ic;
            rpnn = [];
            for kk = 1:length(pname)
                pn = char(get(pname(kk),'string'));
                if ~isempty(pn)
                    present = 0;
                    kkk = 1;
                    while present == 0 & kkk <= NN
                        if ~isempty(findstr(pn,char(der(kkk))))
                            present = 1;
                        end
                        kkk = kkk + 1;
                    end
                    kkk = 1;
                    while present == 0 & kkk <= krex
                        rexval = get(exval(renn(kkk)),'string');
                        if ~isempty(findstr(pn,char(rexval)))
                            present = 1;
                        end
                        kkk = kkk + 1;
                    end
                    if strcmp(flag,'write')
                        while present == 0 & kkk <= NN
                            if ~isempty(findstr(pn,icond{kkk}))
                                present = 1;
                            end
                            kkk = kkk + 1;
                        end
                    end
                    
                    if present
                        rpnn = [rpnn,kk];
                    end
                end
            end
            
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            function writefile(fname, pathname, renn, rpnn)
            
            %% Writes a function M-file for the system.
            
            odesolveset = findobj('name','ODESOLVE Setup');
            sud = get(odesolveset,'user');
            NN = sud.c.NN;
            xname = sud.c.xname;
            tname = sud.c.tname;
            der = sud.c.der;
            exname = sud.h.exname;
            exval = sud.h.exval;
            krex = length(renn);
            rexname = cell(1,krex);
            for kk = 1:krex
                rexname{kk} = get(exname(renn(kk)),'string');
                rexval{kk} = get(exval(renn(kk)),'string');
            end
            pname = sud.h.pname;
            pval = sud.h.pval;
            icond = sud.c.ic;
            krp = length(rpnn);
            for kk = 1:krp
                pn = get(pname(rpnn(kk)),'string');
                rpname{kk} = [',', pn];
                rpns{kk} = pn;
                rpval{kk} = get(pval(rpnn(kk)),'string');
            end  
            
            %% Check if tname is Greek.
            
            if (abs(tname(1)) == 92)
                tname(1) = [];
            end
            
            %% Find a convenient variable name not included in tname, xname,
            %% rpname, or rexname. 
            
            present = 1;
            tee = clock;
            tee = num2str(ceil(10000*tee(6)/6));
            kk=1;
            while present
                if kk == 1
                    uname = 'X';
                else
                    uname = ['X',tee(1,kk-1)]
                end
                present = ~isempty(strmatch(uname,xname));
                if ~present
                    present = ~isempty(strmatch(uname,tname));
                end
                if ~present & krp > 0
                    present = ~isempty(strmatch(uname,rpns));
                end
                if ~present & krex > 0
                    present = ~isempty(strmatch(uname,rexname));
                end
                kk = kk + 1;
            end % while present
            
            %% Start the file.
            
            if krp > 0
                fcnstr = ['function ', uname, 'pr = ',fname,'(',tname,',', uname, rpname{:},')\n\n'];
            else 
                fcnstr = ['function ', uname, 'pr = ',fname,'(',tname,',', uname, ')\n\n'];
            end
            
            cxn{1} = xname{1};
            for kk = 2:NN
                cxn{kk} = [', ', xname{kk}];
            end
            if krp
                cpn{1} = rpns{1};
                for kk = 2:krp
                    cpn{kk} = [', ', rpns{kk}];
                end      
            end
            
            dirstr{1} = '%%Input:\n';
            dirstr{2} = ['%%   Independent variable: ', tname, '\n'];
            dirstr{3} = ['%%   Dependent variable vector: X = [', cxn{:},']\n'];
            if krp == 1
                pdirstr{1} = ['%%   Parameter: ', cpn{:}, '\n'];
                pdirstr{2} = ['%%   ', cpn{:}, ' is a parameter\n'];
            elseif krp > 1
                pdirstr{1} = ['%%   Parameters: ', cpn{:}, '\n'];
                pdirstr{2} = ['%%   ', cpn{:}, ' are parameters\n'];
            end
            outstr = ['%%   Vector of derivatives: ', uname 'pr\n'];
            dirstr{4} = '%%To solve the system, execute\n';
            if krp
                dirstr{5} = ['%%      [t,u] = ode45(@', fname, ', tspan, u0, []', ...
                        rpname{:},')\n'];
            else
                dirstr{5} = ['%%      [t,u] = ode45(@', fname, ', tspan, u0)\n'];
            end
            dirstr{6} = '%%   tspan is the solution interval\n';
            dirstr{7} = '%%   u0 is the initial vector\n';
            commstr = '%%%%%% Created by ODESOLVE\n';
            
            dff = fopen([pathname,fname,'.m'],'w');
            fprintf(dff,fcnstr);
            fprintf(dff,commstr);
            fprintf(dff,'%%\n');
            for kk = 1:3
                fprintf(dff,dirstr{kk});
            end
            if krp
                fprintf(dff,pdirstr{1});
            end
            fprintf(dff,'%%\n');
            fprintf(dff,'%%Output:\n');
            fprintf(dff,outstr);
            fprintf(dff,'%%\n');
            for kk = 4:7
                fprintf(dff,dirstr{kk});
            end
            if krp
                fprintf(dff,pdirstr{2});
            end
            fprintf(dff,'\n\n');
            
            
            
            %% Enter the expressions.
            
            if krex
                for kk = 1:krex
                    exstr = rexval{kk};
                    
                    % Remove all periods
                    
                    l=length(exstr);
                    for ( k = fliplr(findstr('.',exstr)))
                        if (find('*/^' == exstr(k+1)))
                            exstr = [exstr(1:k-1), exstr(k+1:l)];
                        end
                        l=l-1;
                    end
                    
                    % Replace variable names by subscripted uname.  A variable must
                    % be preceded and followed by an operator, unless it is the first
                    % or last symbol in an expression.
                    
                    lopstr = '(+-*/^';
                    ropstr = ')+-*/^';
                    for ll = 1:NN
                        UN = [uname, '(', int2str(ll),')'];
                        XN = char(xname{ll});
                        if length(XN) <= length(exstr)
                            kkk = findstr(exstr,XN);
                            Le = length(XN);
                            while ~isempty(kkk)
                                Lde = length(exstr);
                                TP = 1; % True presence at this location is assumed.
                                kkkk = max(kkk);
                                if kkkk + Le < Lde & isempty(find(ropstr == exstr(kkkk+Le)))
                                    TP = 0;
                                end
                                if kkkk > 1 & isempty(find(lopstr == exstr(kkkk-1)))
                                    TP = 0;
                                end
                                if TP
                                    exstr = [exstr(1:kkkk-1),UN,exstr(kkkk+Le:Lde)];
                                end
                                if Le < kkkk
                                    kkk = findstr(exstr(1:kkkk-1),XN);
                                else
                                    kkk = [];
                                end
                            end
                        end
                    end
                    eqnstr = [rexname{kk}, ' = ', char(exstr),';\n'];
                    fprintf(dff,eqnstr);
                end
                fprintf(dff,'\n');
            end     
            
            %% Enter the derivatives.
            
            for kk = 1:NN
                derstr = char(der{kk});
                derstr(find(abs(derstr)==32))=[];   % Remove blanks.
                
                % Remove all periods
                
                l=length(derstr);
                for ( k = fliplr(findstr('.',derstr)))
                    if (find('*/^' == derstr(k+1)))
                        derstr = [derstr(1:k-1), derstr(k+1:l)];
                    end
                    l=l-1;
                end
                
                
                % Replace variable names by subscripted uname.  A variable must
                % be preceded and followed by an operator, unless it is the first
                % or last symbol in an expression.
                
                lopstr = '(+-*/^';
                ropstr = ')+-*/^';
                for ll = 1:NN
                    UN = [uname, '(', int2str(ll),')'];
                    XN = char(xname{ll});
                    if length(XN) <= length(derstr)
                        kkk = findstr(derstr,XN);
                        Le = length(XN);
                        while ~isempty(kkk)
                            Lde = length(derstr);
                            TP = 1; % True presence at this location is assumed.
                            kkkk = max(kkk);
                            if kkkk + Le < Lde & isempty(find(ropstr == derstr(kkkk+Le)))
                                TP = 0;
                            end
                            if kkkk > 1 & isempty(find(lopstr == derstr(kkkk-1)))
                                TP = 0;
                            end
                            if TP
                                derstr = [derstr(1:kkkk-1),UN,derstr(kkkk+Le:Lde)];
                            end
                            if Le < kkkk
                                kkk = findstr(derstr(1:kkkk-1),XN);
                            else
                                kkk = [];
                            end
                        end
                    end
                end
                eqnstr = [derstr,';'];
                if NN == 1
                    eqnstr = [uname, 'pr = ', eqnstr];
                else
                    if kk == 1
                        eqnstr = [uname, 'pr = [', eqnstr];
                    else
                        eqnstr = [blanks(7), eqnstr];
                    end
                    if kk < NN
                        eqnstr = [eqnstr, '...\n'];
                    else
                        eqnstr = [eqnstr, '];\n'];
                    end
                end
                fprintf(dff,eqnstr);
            end   
            fclose(dff);