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

% DFIELD6       is an interactive tool for studying single first order
%       differential equations.  When DFIELD6 is executed, a DFIELD6 Setup
%       window is opened.  The user may enter the differential
%       equation and specify a display window using the interactive
%       controls in the Setup window.
%
%       When the Proceed button is pressed on the Setup window, the DF
%       Display window is opened.  At first this window displays a
%       direction line field for the differential equation.  When the
%       mouse button is depressed in the DFIELD6 Display window, the
%       solution to the differential equation with that initial
%       condition is calculated and plotted.
%
%       Other options are available in the Options menu.  These are
%       fairly self explanatory.  The Settings option allows the user
%       to change several parameters.  Included here are the
%       possibilities of using a vector field instead of the default
%       line field, and of changing the number of field points computed
%       and displayed.

%  Copywright (c) 1995, 1997, 1998, 1999, 2000 
%                 John C. Polking Rice University
%  Last modified: January 21, 2001


startstr = 'dfield6';

if nargin < 1
  action ='initialize';
end


% disp(action)


if strcmp(action,'initialize')

  % First we make sure that there is no other copy of DFIELD6
  % running, since this causes problems.        
  
  dfh = findobj('tag','DFIELD6');
  if ~isempty(dfh);
    qstring = {'There are some open DFIELD6 figures (perhaps invisible).  ';...
    'What do you want to do?'};
    tstring = 'Only one copy of DFIELD6 can be open at one time.';
    b1str = 'Restart DFIELD6.';
    b2str = 'Just close those figures.';
    b3str = 'Do nothing.';
    answer = questdlg(qstring,tstring,b1str,b2str,b3str,b1str);
    if strcmp(answer,b1str)
      delete(dfh);
      eval(startstr);return
    elseif strcmp(answer,b2str)
      delete(dfh);return
    else
      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 DFIELD6.
  
  oldfiles = dir('dftp*.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 DFIELD6')
      delete(fn)
    else
      kk = [kk;k];
    end
  end
  oldfiles = dir([tempdir,'dftp*.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 DFIELD6')
      delete(fn)
    else
      kk = [kk;k];
    end
  end
  if ~isempty(kk)
    fprintf('The following files\n');
    for j = 1:length(kk)
      fprintf(['     ',oldfiles(kk(j)).name,'\n']);
    end
    fprintf('may have been created by an older version of DFIELD.\n');
    fprintf('If so they should be deleted.\n\n');
  end
  
  style = 'white';
  ssize = get(0,'defaultaxesfontsize');
  npts = 20;
  solver = 'Dormand-Prince';
  if exist('dfstart','file')
    H = dfstart;
    if ~isempty(H)
      if isfield(H,'style')
        style = H.style;
      end
      if isfield(H,'size')
        ssize = H.size;
      end
      if isfield(H,'npts')
        npts = H.npts;
      end
      if isfield(H,'solver')
        solver = H.solver;  
      end
    end
  end
if get(0,'screendepth') == 1
  style = 'bw';
end
ud.ssize = ssize;
  comp = computer;
  if strcmp(comp,'PCWIN')
    ud.fontsize = 0.8*ud.ssize;
  else
    ud.fontsize = ud.ssize;
  end
  system.xname = 'x';
  system.tname = 't';
  system.der = ' x^2 - t';
  system.wind = [-2 10 -4 4];
  system.pname1 = '';
  system.pname2 = '';
  system.pval1 = '';
  system.pval2 = '';
  
  ud.c = system;        % Changed values.
  ud.o = system;        % Original values.
  ud.h = system;        % This will be the handles in the
                        % setup window.
  ud.fieldtype = 'lines';
  ud.npts = npts;  
  ud.style = style;
  ud.settings.magn = 1.25;
  ud.settings.refine = 4;
  ud.settings.tol = 1e-6;
  ud.settings.solver = solver;
  ud.settings.stepsize = 0.1;
  ud.settings.hmax = 0;
  ud.settings.speed = 100;
  ud.systems = system;
  switch style
  case 'black'
    color.temp = [1 0 0]; % red for temporary orbits
    color.orb = [1 1 0];  % yellow for orbits
    color.arrows = .5*[1 1 1];  % gray for arrows
    color.pline = [1 1 1];    
    color.level = [1,.5,.5];
   
  case 'white'
    color.temp = [1 0 0]; % red for temporary orbits
    color.orb = [0 0 1];  % blue for orbits
    color.arrows = .7*[1 1 1]; % gray for arrows
    color.pline = [0 0 0];    
    color.level = 0.8*[.9,.5,.8];
    
  case 'test'
    color.temp = [1 0 0]; % red for temporary orbits
    color.orb = [0 0 1];  % blue for orbits
    color.arrows = .7*[1 1 1]; % gray for arrows
    color.pline = [0 0 0];    
    color.level = [1,.5,.5];
    
  case 'display'
    color.temp = [1 0 0]; % red for temporary orbits
    color.orb = [0 0 1];  % blue for orbits
    color.arrows = .4*[1 1 1]; % gray for arrows
    color.pline = [0 0 0];   
    color.level = [255 150 10]/255;%  0.8*[.9,.5,.8];
    
  case 'bw'
    color.temp = [1 1 1]; % white for everything
    color.orb = [1 1 1];  
    color.arrows = [1 1 1];
    color.pline = [1 1 1];    
    color.level = [1,1,1];
  end
  ud.color = color;
  
  dfset = figure('name','DFIELD6 Setup','numb','off',...
      'tag','DFIELD6','visible','off',...
      'user',ud);
  
  dfield6('figdefault',dfset);
  frame(1) = uicontrol('style','frame','visible','off');
  eq(1)=uicontrol('style','text',...
      'horizon','center',...
      'string','The differential equation.','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 = 0;             % Left margin of the frames.
  frsep = 0;            % Separation between frames.
  separation = texth;   % Separation between boxes.

  disfrw = varw+equalw+eqlength+10;
  
  dfigwidth =2*left + disfrw;   % Width of the figure.
  dfigurebot = 60;      % Bottom of the figure.
  
  buttw = dfigwidth/3;
  qwind = [0,frsep,buttw,separation];   % Quit button
  rwind = [buttw,frsep,buttw,separation];       % Revert "
  pwind = [2*buttw,frsep,buttw,separation];     % Proceed "
  
  disfrbot = 2*frsep + separation;      % Display frame.
  disfrht = 3*separation + 10;
  disfrwind = [left, disfrbot, disfrw, disfrht];
    
  pfrbot = disfrbot + disfrht + frsep;
  pfrw = disfrw;
  pfrht = separation + 10;
  pfrwind = [left, pfrbot, pfrw, pfrht];
  
  defrbot = pfrbot + pfrht + frsep;     % Equation frame.
  defrw = pfrw;
  defrht = 3*separation + 15;
  defrwind = [left, defrbot, defrw, defrht];
    
  dfigureheight = defrbot + defrht +frsep;      % Height of the figure.
  
  set(dfset,'pos',[30 dfigurebot dfigwidth dfigureheight]);
  
  set(frame(1),'pos',defrwind);
  
  
  xname=[
    'ud = get(gcf,''user'');'...
        'Xname=get(ud.h.xname,''string'');'...
        'set(ud.h.twind(3),''string'','...
           '[''The minimum value of '',Xname,'' = '']);'...
        'set(ud.h.twind(4),''string'','...
           '[''The maximum value of '',Xname,'' = '']);'...
        'ud.c.xname = Xname;'...
        'ud.flag = 0;'...
        'set(gcf,''user'',ud);'];
  
  tname=[
    'ud = get(gcf,''user'');'...
        'Tname=get(ud.h.tname,''string'');'...
        'set(ud.h.twind(1),''string'','...
          '[''The minimum value of '',Tname,'' = '']);'...
        'set(ud.h.twind(2),''string'','...
          '[''The maximum value of ''Tname,'' = '']);'...
        'ud.c.tname = Tname;'...
        'ud.flag = 0;'...
        'set(gcf,''user'',ud);'];
  
  der =[
    'ud = get(gcf,''user'');'...
        'ud.c.der = get(ud.h.der,''string'');'...
        'ud.flag = 0;'...
        'set(gcf,''user'',ud);'];
    
  equationbot = defrbot + 5;
  eqlabelbot = equationbot + 2*separation;
  derbot = equationbot + separation + 5;                % Bottom of equation.
  tbot = equationbot;                           % Bottom of ind. var.
  
  lablen =200*rr;
  
  eqlableft = (dfigwidth-lablen)/2;
  eqleft = left + 5;
  fudge = 0.15*separation;
  set(eq(1),'pos',[eqlableft eqlabelbot lablen texth]);
  tcolor = get(gcf,'defaultuicontrolbackgroundcolor');
  %if strcmp(style,'test')
  %   tcolor = [1 1 1];
  %end
  
  ud.h.xname=uicontrol('pos',[eqleft, derbot, varw, texth],...
     'style','edit',...
     'horizon','right',...
     'string',ud.o.xname,...
     'call',xname,...
     'visible','off',...
     'backgroundcolor',tcolor);
  
  eq(2) = uicontrol('style','text',...
      'pos',[eqleft+varw derbot-fudge equalw texth],...
      'horizon','center','string',''' = ','visible','off',...
      'backgroundcolor',tcolor);
  
  ud.h.der=uicontrol('pos',[eqleft+varw+equalw derbot eqlength texth],...
      'string',ud.o.der,...
      'horizon','left','style','edit',...
      'call',der,'visible','off','backgroundcolor',tcolor);
  
  eq(3) = uicontrol('style','text',...
      'horizon','right','string','The independent variable is ',...
      'visible','off','backgroundcolor',tcolor);
 
  ext = get(eq(3),'extent');
  tnw = ext(3)+10;
  tvarl = eqleft + tnw;
  set(eq(3),'pos',[eqleft,tbot-fudge,tnw,texth]);
  
  ud.h.tname = uicontrol('pos',[tvarl tbot varw texth],'style','edit',...
     'horizon','left','string',ud.o.tname,'call',tname,...
     'visible','off','backgroundcolor',tcolor);
  
  pframe = uicontrol('style','frame','pos',pfrwind,'visible','on');
  
  p1 = [
    'ud = get(gcf,''user'');'...
        'ud.c.pname1=get(ud.h.pname1,''string'');'...
        'ud.c.pval1=get(ud.h.pval1,''string'');'...
        'ud.flag = 0;'...
        'ud.c.name = '''';',...
        'set(gcf,''user'',ud);'];
  
  p2 = [
    'ud = get(gcf,''user'');'...
        'ud.c.pname2=get(ud.h.pname2,''string'');'...
        'ud.c.pval2=get(ud.h.pval2,''string'');'...
        'ud.flag = 0;'...
        'ud.c.name = '''';',...
        'set(gcf,''user'',ud);'];

  pnamew = 30*rr;
  pvalw = 55*rr;
  peqw = 8*rr;
  pbot = pfrbot + 5;

  
  paratit=uicontrol('style','text',...
      'horizon','center',...
      'string','Parameters/expressions:',...
      'visible','off','backgroundcolor',tcolor);
  
  ext = get(paratit,'extent');
  paratitw = ext(3);
  pos = [eqleft pbot-fudge paratitw texth];
  set(paratit,'pos',pos);
  pleft1 = eqleft + paratitw + 3;
  peqleft1 = pleft1 + pnamew;
  pvleft1 = peqleft1 + peqw;
  pleft2 = dfigwidth - 10 - pnamew - pvalw - peqw;
  peqleft2 = pleft2 + pnamew;
  pvleft2 = peqleft2 + peqw;    

  ud.h.pname1 = uicontrol('style','edit',...
      'pos',[pleft1 pbot pnamew texth],...
      'horizon','right','string',ud.o.pname1,...
      'call',p1,'visible','off',...
      'backgroundcolor',tcolor);
  
  equal(1) = uicontrol('style','text',...
      'pos',[peqleft1 pbot-fudge peqw texth],...
      'horizon','center','string','=','visible','off',...
      'backgroundcolor',tcolor);
  
  ud.h.pval1 = uicontrol('style','edit',...
      'pos',[pvleft1 pbot pvalw texth],...
      'string',ud.o.pval1,'call',p1,'visible','off',...
      'backgroundcolor',tcolor);
  
  ud.h.pname2 = uicontrol('style','edit',...
      'pos',[pleft2 pbot pnamew texth],...
      'horizon','right','string',ud.o.pname2,...
      'call',p2,'visible','off',...
      'backgroundcolor',tcolor);
  
  equal(2) = uicontrol('style','text',...
      'pos',[peqleft2 pbot-fudge peqw texth],...
      'horizon','center','string','=','visible','off',...
      'backgroundcolor',tcolor);
  
  ud.h.pval2 = uicontrol('style','edit',...
      'pos',[pvleft2 pbot pvalw texth],...
      'string',ud.o.pval2,'call',p2,'visible','off',...
      'backgroundcolor',tcolor);
  
  
  frame(2) = uicontrol('style','frame','pos',disfrwind,'visible','off');
  
  w1 = [
    'ud = get(gcf,''user'');'...
        'nnn =  str2num(get(ud.h.wind(1),''string''));'...
        'if isempty(nnn),',...
        '  set(ud.h.wind(1),''string'',''?'');',...
        '  nnn = NaN;',...
        'end,',...
        'ud.c.wind(1) = nnn;',...           
        'set(gcf,''user'',ud);'];       
  
  w2 = [
    'ud = get(gcf,''user'');'...
        'nnn =  str2num(get(ud.h.wind(2),''string''));'...
        'if isempty(nnn),',...
        '  set(ud.h.wind(2),''string'',''?'');',...
        '  nnn = NaN;',...
        'end,',...
        'ud.c.wind(2) = nnn;',...           
        'set(gcf,''user'',ud);'];       
  
  w3 = [
    'ud = get(gcf,''user'');'...
        'nnn =  str2num(get(ud.h.wind(3),''string''));'...
        'if isempty(nnn),',...
        '  set(ud.h.wind(3),''string'',''?'');',...
        '  nnn = NaN;',...
        'end,',...
        'ud.c.wind(3) = nnn;',...           
        'set(gcf,''user'',ud);'];       
  
  w4 = [
    'ud = get(gcf,''user'');'...
        'nnn =  str2num(get(ud.h.wind(4),''string''));'...
        'if isempty(nnn),',...
        '  set(ud.h.wind(4),''string'',''?'');',...
        '  nnn = NaN;',...
        'end,',...
        'ud.c.wind(4) = nnn;',...           
        'set(gcf,''user'',ud);'];       
  
  
  winbot1 = disfrbot + disfrht - 5 - separation;
  winbot2 = winbot1 - separation;
  winbot3 = winbot2 - separation;
  
  windw = 40*rr;  
  twindw = (disfrw - 10)/2-windw;
  twindl = eqleft + twindw + windw;

  dwind = uicontrol('style','text',...
      'pos',[eqleft winbot1 disfrw-10 texth],...
      'horizon','center',...
      'string','The display window.','visible','off',...
      'backgroundcolor',tcolor);
  
  % ud.h.twind contains the handles to the text windows, and ud.h.wind
  % contains the handles to the edit windows.
  
  ud.h.twind(1) = uicontrol('style','text',...
      'pos',[eqleft winbot2-fudge twindw texth],...
      'horizon','right',...
      'string',['The minimum value of ',ud.o.tname,' = '],...
      'visible','off','backgroundcolor',tcolor);
  
  ud.h.wind(1) = uicontrol('style','edit',...
      'pos',[eqleft+twindw winbot2 windw texth],...
      'string',num2str(ud.o.wind(1)),...
      'call',w1,'visible','off',...
      'backgroundcolor',tcolor);
  
  ud.h.twind(2) = uicontrol('style','text',...
      'pos',[eqleft winbot3-fudge twindw texth],...
      'horizon','right',...
      'string',['The maximum value of ',ud.o.tname,' = '],...
      'visible','off','backgroundcolor',tcolor);
  
  ud.h.wind(2) = uicontrol('style','edit',...
      'pos',[eqleft+twindw winbot3 windw texth],...
      'string',num2str(ud.o.wind(2)),...
      'call',w2,'visible','off',...
      'backgroundcolor',tcolor);
  
  ud.h.twind(3)= uicontrol('style','text',...
      'pos',[twindl winbot2-fudge twindw texth],...
      'horizon','right',...
      'string',['The minimum value of ',ud.o.xname,' = '],...
      'visible','off','backgroundcolor',tcolor);
  
  ud.h.wind(3) = uicontrol('style','edit',...
      'pos',[twindl+twindw winbot2 windw texth],...
      'string',num2str(ud.o.wind(3)),...
      'call',w3,'visible','off',...
      'backgroundcolor',tcolor);
  
  ud.h.twind(4) = uicontrol('style','text',...
      'pos',[twindl winbot3-fudge twindw texth],...
      'horizon','right',...
      'string',['The maximum value of ',ud.o.xname,' = '],...
      'visible','off','backgroundcolor',tcolor);
  
  ud.h.wind(4) = uicontrol('style','edit',...
      'pos',[twindl+twindw winbot3 windw texth],...
      'string',num2str(ud.o.wind(4)),...
      'call',w4,'visible','off',...
      'backgroundcolor',tcolor);
  
  
  
  butt(1) = uicontrol('style','push',...
      'pos',qwind,...
      'string','Quit','call','dfield6(''quit'')','visible','off');
  
  butt(2) = uicontrol('style','push',...
      'pos',rwind,...
      'string','Revert','call','dfield6(''revert'')','visible','off');
  
  butt(3) = uicontrol('style','push',...
      'pos',pwind,...
      'string','Proceed','call','dfield6(''proceed'')','visible','off');

  hhsetup = get(0,'showhiddenhandles');
  set(0,'showhiddenhandles','on');
  
  mefile = findobj(dfset,'label','&File');
  meedit = findobj(gcf,'label','&Edit');
  metools = findobj(gcf,'label','&Tools');
  meview = findobj(gcf,'label','&View');
  meinsert = findobj(gcf,'label','&Insert');
  delete([metools,meview,meinsert]);
  
  % 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');
  meproceed = uimenu(mefile,'label','Proceed','call',...
      'dfield6(''proceed'')','separator','off',...
      'accelerator','G','pos',1);
  merevert = uimenu(mefile,'label','Revert','call',...
      'dfield6(''revert'')','separator','off','pos',2);
  set(mepset,'vis','on','pos',3);
  set(meprev,'vis','on','label','Page Pre&view...','pos',4);
  set(meexp,'vis','on','pos',5,'separator','off');
  merestart = uimenu(mefile,'label','Restart DFIELD6','call',...
      'dfield6(''restart'')','separator','on','pos',6);
  mequit = uimenu(mefile,'label','Quit DFIELD6','call',...
      'dfield6(''quit'')','separator','on','pos',7);
  
  set(get(meedit,'child'),'vis','off');
  meclrf = uimenu(meedit,'label','Clear equations',...
      'call',['ud = get(gcf,''user'');h = ud.h;',...
        'set([h.xname,h.der,h.tname,h.der],''string'','''');'],...
      'accelerator','E');
  meclrp = uimenu(meedit,'label','Clear parameters',...
      'call',['ud = get(gcf,''user'');h = ud.h;',...
        'set([h.pname1,h.pname2,h.pval1,h.pval2],''string'','''');'],...
      'accelerator','N');
  meclrwind = uimenu(meedit,'label','Clear display window',...
      'call',['ud = get(gcf,''user'');',...
        'set(ud.h.wind,''string'','''');'],...
      'accelerator','D');
  meclrall = uimenu(meedit,'label','Clear all',...
      'call',['ud = get(gcf,''user'');h = ud.h;',...
        'set(h.wind,''string'','''');',...
        'set([h.xname,h.der,h.tname,h.der,h.pname1,h.pname2,h.pval1,h.pval2],''string'','''');'],...
      'accelerator','A',...
      'separator','on');
  set(0,'showhiddenhandles',hhsetup);
  ud.flag = 0;
  
  % Record the handles in the User Data of the Set Up figure.
  
  set(dfset,'user',ud);
  hhhh = findobj(dfset,'type','uicontrol');
  set(hhhh,'units','normal')
  
  set(dfset,'visible','on','resize','on');
  set(get(dfset,'children'),'visible','on');
  
  
elseif strcmp(action,'revert')
  
  ud = get(gcf,'user');
  ud.c = ud.o;
  syst = ud.o;
  xname = syst.xname;
  tname = syst.tname;
  set(ud.h.xname,'string',xname);
  set(ud.h.tname,'string',tname);
  set(ud.h.der,'string',syst.der);
  set(ud.h.twind(1),'string',['The minimum value of ',tname,' = ']);
  set(ud.h.twind(2),'string',['The maximum value of ',tname,' = ']);
  set(ud.h.twind(3),'string',['The minimum value of ',xname,' = ']);
  set(ud.h.twind(4),'string',['The maximum value of ',xname,' = ']);
  set(ud.h.wind(1),'string',num2str(syst.wind(1)));
  set(ud.h.wind(2),'string',num2str(syst.wind(2)));
  set(ud.h.wind(3),'string',num2str(syst.wind(3)));
  set(ud.h.wind(4),'string',num2str(syst.wind(4)));
  set(ud.h.pname1,'string',syst.pname1);
  set(ud.h.pname2,'string',syst.pname2);
  set(ud.h.pval1,'string',syst.pval1);
  set(ud.h.pval2,'string',syst.pval2);
   
  set(gcf,'user',ud);
  
elseif strcmp(action,'proceed')

  % Proceed connects Setup with the Display window.

  dfset = gcf;
  sud = get(dfset,'user');
  sud.o = sud.c;
  set(dfset,'user',sud);

  % Some error checking that has to be done no matter what.
  
  WINvect = sud.c.wind;
  if any(isnan(WINvect))
    errmsg = ['One of the entries defining the display window ',...
          'is not a number.'];
    fprintf('\a')
    errordlg(errmsg,'DFIELD6 error','on');
    return
  end
  xstr = sud.c.xname;
  if isempty(xstr)
    errmsg = 'The dependent variable needs a name.';
    fprintf('\a')
    errordlg(errmsg,'DFIELD6 error','on');
    return
  end
  tstr = sud.c.tname;
  if isempty(tstr)
    errmsg = 'The independent variable needs a name.';
    fprintf('\a')
    errordlg(errmsg,'DFIELD6 error','on');
    return
  end
  if WINvect(2)<= WINvect(1)
    errmsg = ['The minimum value of ', tstr,...
          ' must be smaller than the maximum value.'];
    fprintf('\a')
    errordlg(errmsg,'DFIELD6 error','on');
    return
  end
  if WINvect(4)<= WINvect(3)
    errmsg = ['The minimum value of ', xstr,...
          ' must be smaller than the maximum value.'];
    fprintf('\a')
    errordlg(errmsg,'DFIELD6 error','on');
    return
  end

  %  sud.flag = 0 if this is the first time through for this equation, 
  %  sud.flag = 1 if only the window dimensions have been changed.
  
  % If sud.flag == 1 we only have to update things.
  
  if (sud.flag == 1)
    dfdisp = findobj('name','DFIELD6 Display');
    dud = get(dfdisp,'user');
    aud = get(dud.axes,'user');
    tstr = get(get(dud.axes,'title'),'string');
    wind = sud.c.wind(:);
    hmax = dud.settings.hmax;
    if (~all(wind == dud.syst.wind(:)))
      dwind = [wind(1); wind(3); -wind(2); -wind(4)];
      DY = [wind(2)-wind(1); wind(4)-wind(3)];
      hmax = min(hmax,DY(1)/4);
      aud.DY = DY;
      aud.cwind = wind - dud.settings.magn*[DY(1);-DY(1);DY(2);-DY(2)];
      set(dud.axes,'user',aud);
    end
    dud.syst = sud.c;
    dud.settings = sud.settings;
    dud.settings.hmax = hmax;
    set(dfdisp,'user',dud);
    dfield6('dirfield',dfdisp);
    
  else
    sud.flag = 1;
    set(dfset,'user',sud);      
    Arrflag = sud.fieldtype;
    NumbFPts = sud.npts;
    Xname = sud.c.xname;
    Tname = sud.c.tname;
    derivstr = sud.c.der;
  
    pname1 = sud.c.pname1;
    pname2 = sud.c.pname2;
    parav1 = sud.c.pval1;
    parav2 = sud.c.pval2;
    
    % First remove the blanks.
    
    derivstr(find(abs(derivstr)==32))=[];
    if ~isempty(parav1)
      parav1(find(abs(parav1)==32))=[];
    end
    if ~isempty(parav2)
      parav2(find(abs(parav2)==32))=[];
    end
    
    % Next remove the periods inserted by users attempting to make the 
    % function array smart.
    
    l=length(derivstr);
    for ( k = fliplr(findstr('.',derivstr)))
      if (find('*/^' == derivstr(k+1)))
        derivstr = [derivstr(1:k-1), derivstr(k+1:l)];
      end
      l=l-1;
    end
    l=length(parav1);
    for ( k = fliplr(findstr('.',parav1)))
      if (find('*/^' == parav1(k+1)))
        parav1 = [parav1(1:k-1), parav1(k+1:l)];
      end
      l=l-1;
    end
    l=length(parav2);
    for ( k = fliplr(findstr('.',parav2)))
      if (find('*/^' == parav2(k+1)))
        parav2 = [parav2(1:k-1), parav2(k+1:l)];
      end
      l=l-1;
    end
        
    % Build strings for the title.
    txderstr = derivstr;
  
    kxder = find(abs(txderstr)==42);  % Get rid of *s
    txderstr(kxder)=' '*ones(size(kxder));
    txderstr = strrep(txderstr,'-',' - ');  % Extra spaces
    txderstr = strrep(txderstr,'+',' + ');
    if (abs(txderstr(1)) == 32)  % Get rid of starting space
      txderstr = txderstr(2:length(txderstr));
    end
    tstr = [Xname,' '' = ', txderstr];
    
    pstring = cell(2,1);
    pstring{1} = '';
    pstring{2} = '';
    
    if ~isempty(parav1)
      tp1str = parav1;
      
      kxder = find(abs(tp1str)==42);  % Get rid of *s
      tp1str(kxder)=' '*ones(size(kxder));
      tp1str = strrep(tp1str,'-',' - ');  % Extra spaces
      tp1str = strrep(tp1str,'+',' + ');
      if (abs(tp1str(1)) == 32)  % Get rid of starting space
        tp1str = tp1str(2:length(tp1str));
      end
      pstring{1} = [pname1,' = ', tp1str];
    end
    
    if ~isempty(parav2)
      tp2str = parav2;
      
      kxder = find(abs(tp2str)==42);  % Get rid of *s
      tp2str(kxder)=' '*ones(size(kxder));
      tp2str = strrep(tp2str,'-',' - ');  % Extra spaces
      tp2str = strrep(tp2str,'+',' + ');
      if (abs(tp2str(1)) == 32)  % Get rid of starting space
        tp2str = tp2str(2:length(tp2str));
      end
      pstring{2} = [pname2,' = ', tp2str];
    end
    
    % Get ready to do some error trapping.
    
    SS = warning;
    warning off
    XxXxXx = WINvect(1) + rand*(WINvect(2)-WINvect(1));
    TtTtTt = WINvect(3) + rand*(WINvect(4)-WINvect(3));
    err = 0;
    
    % Now we remove the backslashes (\) used to get Greek into the
    % labels. 
    txname = Xname;
    ttname = Tname;   
    derivstr(find(abs(derivstr)==92))=[];
    Xname(find(abs(Xname)==92))=[];
    Tname(find(abs(Tname)==92))=[];
    
    eval([Xname,'=XxXxXx;'],'err = 1;');
    if err
      errmsg = ['"',xstr, '" is not a valid variable name in MATLAB.'];
      fprintf('\a')
      errordlg(errmsg,'DFIELD6 error','on');
      warning(SS)
      return
    end
    err = 0;
    eval([Tname,'=TtTtTt;'],'err = 1;');
    if err
      errmsg = ['"',tstr, '" is not a valid variable name in MATLAB.'];
      fprintf('\a')
      errordlg(errmsg,'DFIELD6 error','on');
      warning(SS)
      return
    end
    
    
    % Replace the parameters/expressions with their values.
    
    pflag = zeros(1,2);
    perr = [];
    if ~isempty(pname1)
      if isempty(parav1) 
        perr = sud.h.pval1;
      else 
        if isempty(str2num(parav1)) % This is an expression.
          tparav1 = parav1;
          tparav1(find(abs(tparav1)==92))=[];
          err = 0; pval = '';
          eval(['pval = ',tparav1,';'],'err=1;');pval;
          if (err | isempty(pval))
            errmsg = ['The expression for ',pname1,' is not valid.'];
            fprintf('\a')
            errordlg(errmsg,'DFIELD6 error','on');
            warning(SS)
            return
          end
        end
        xderivstr = dfield6('paraeval',pname1,parav1,derivstr);
        if ~strcmp(xderivstr,derivstr)
          pflag(1) = 1;
          derivstr = xderivstr;
        end
      end
    end

    if ~isempty(pname2)
      if isempty(parav2) 
        perr = sud.h.pval2;
      else 
        if isempty(str2num(parav2)) % This is an expression.
          tparav2 = parav2;
          tparav2(find(abs(tparav2)==92))=[];
          err = 0; pval = '';
          eval(['pval = ',tparav2,';'],'err=1;');pval;
          if (err | isempty(pval))
            errmsg = ['The expression for ',pname2,' is not valid.'];
            fprintf('\a')
            errordlg(errmsg,'DFIELD6 error','on');
            warning(SS)
            return
          end
        end
        xderivstr = dfield6('paraeval',pname2,parav2,derivstr);
        if ~strcmp(xderivstr,derivstr)
          pflag(2) = 1;
          derivstr = xderivstr;
        end
      end
    end

    % We have to make the derivative strings array smart.
    
    l = length(derivstr);
    for (k=fliplr(find((derivstr=='^')|(derivstr=='*')|(derivstr=='/'))))
      derivstr = [derivstr(1:k-1) '.' derivstr(k:l)];
      l = l+1;
    end
    
    % Some more error trapping.
    
    eval(['res = ',derivstr, ';'],'err = 1;');
    if err | isempty(res)
       if isempty(perr)
          errmsg = ['The differential equation ',...
                'is not entered correctly.'];
       else
          errstr1 = 'The differential equation does not evaluate correctly.';
          errstr2 = 'At least one of the parameter values is not a number.';
          errmsg = str2mat(errstr1,errstr2);
          set(perr,'string','?');
       end
       fprintf('\a')
       errordlg(errmsg,'DFIELD6 error','on');
       dfield6('restart');
       warning(SS);
       return
    end
    warning(SS);
    
    % Build a new derivative m-file.
    
    tee = clock;
    tee = ceil(tee(6)*100);
    dfcn=['dftp',num2str(tee)];
    fcnstr = ['function YYyYypr = ',dfcn,'(TtTt,YyYy)\n\n'];
    commstr = '%%%% Created by DFIELD6\n\n';
    varstr = [Xname,' = YyYy;', Tname,' = TtTt;\n\n'];
    lenstr = ['l = length(YyYy);\n'];
    derstr1 = ['YYyYypr = ', derivstr,';\n'];
    derstr2 = ['if (length(YYyYypr) < l) YYyYypr = YYyYypr*ones(1,l);end\n'];

    dff = fopen([tempdir,dfcn,'.m'],'w');
    fprintf(dff,fcnstr);
    fprintf(dff,commstr);
    fprintf(dff,varstr);
    fprintf(dff,lenstr);
    fprintf(dff,derstr1);
    fprintf(dff,derstr2);
    fclose(dff);


    % Find DFIELD6 Display if it exists.        
    % If DFIELD6 Display exists, update it.  If it does not build it.
    
    dfdisp = findobj('name','DFIELD6 Display');
    if (~isempty(dfdisp))
      figure(dfdisp);
      dud = get(dfdisp,'user');
      dud.syst = sud.c;
      dud.settings = sud.settings;
      dfcnn = dud.function;
      if exist(dfcnn) == 2, delete([dfcnn,'.m']);end
      
    else
      dfdisp = figure('name','DFIELD6 Display',...
          'numb','off','interrupt','on',...
          'visible','off','tag','DFIELD6');
      dfield6('figdefault',dfdisp);
      dud = get(dfdisp,'user'); 
      dud.syst = sud.c;
      dud.settings = sud.settings;
      fs = dud.fontsize;
      ssize = dud.ssize;
      r = ssize/10;
      dfaxw = 437*1.2;  % Default axes width
      dfaxh = 315*1.2;  % Default axes height
      dfaxl = 45*1.2;   % Default axes left
      buttw = 40*1.2;     % Default button width
      uni = get(0,'units');
      set(0,'units','pixels');
      ss = get(0,'screensize');
      style = sud.style;
      set(0,'units',uni);
      
      nframeh = 70;   % Default notice frame height
      dfaxb = 4+nframeh+35;

      titleh = 30;
      bottomedge = 38;
      dfdh = dfaxh + nframeh + titleh + bottomedge;

      sw = ss(3);sh = ss(4);
      bottom = 20;
      
      
      if r*dfdh > sh - bottom -35;
         r = (sh-bottom-35)/dfdh;
         fs = 10*r;
         lw = 0.5*r;
         set(gcf,'defaultaxesfontsize',fs,'defaultaxeslinewidth',lw);
         set(gcf,'defaulttextfontsize',fs);
         set(gcf,'defaultlinelinewidth',lw);
         set(gcf,'defaultuicontrolfontsize',fs*0.9);
      end
      
      % Set up the bulletin window.
      
      nframe = uicontrol('style','frame','visible','on');
      nstr = {'More';'than';'five';'lines';'of text'};
      dud.notice = uicontrol('style','text',...
          'horiz','left',...
          'string',nstr,'visible','on');
      ext = get(dud.notice,'extent');
      nframeh = ext(4)+2;
      
      titleh = r*titleh;
      dfaxl = r*dfaxl;
      buttw = r*buttw;
      butth = fs+20*r;
 
      dfaxw = r*dfaxw;
      dfaxh = r*dfaxh;
      dfaxb = nframeh+r*bottomedge;
      buttl = dfaxl + dfaxw + 5;
      buttsep = (dfaxh - butth)/2;
      
      % Set up the plot axes.
    
      dud.axes = axes('units','pix',...
         'position',[dfaxl,dfaxb,dfaxw,dfaxh],...
         'next','add',...
         'box','on',...
         'interrupt','on',...
         'xgrid','on',...
         'ygrid','on',...
         'drawmode','fast',...
         'visible','off',...
         'tag','display axes');
      
      % Finish the positions.
      
      dfdw = buttl + buttw +5;
      dfdh = dfaxb+dfaxh+titleh;
      set(nframe,'pos',[10,1,dfdw-20,nframeh]);
      set(dud.notice,'pos',[15,1,dfdw-30,nframeh-2],...
          'string',{' ';' ';' ';' ';' '});
      
      dfdleft = max((sw-dfdw)/2,sw-dfdw-60);
      dfdbot = sh - dfdh - 70;
      dfdpos = [dfdleft,dfdbot,dfdw,dfdh];
      set(dfdisp,'resize','on');
      set(dfdisp,'pos',dfdpos);
      
      Arrflag = sud.fieldtype;
      
      % Set up the buttons and the menu.
      
      stopstr = 'aud = get(gca,''user'');aud.stop = 4;set(gca,''user'',aud);';
      
      dbutt(1) = uicontrol('style','push',...
          'pos',[buttl,dfaxb+2*buttsep,buttw,butth],...
          'string','Stop','call',stopstr,'vis','off','tag','stop');
      
      dbutt(2) = uicontrol('style','push',...
          'pos',[buttl,dfaxb,buttw,butth],...
          'string','Quit',...
          'call','dfield6(''quit'')','visible','off');
      
      dbutt(3) = uicontrol('style','push',...
          'pos',[buttl,dfaxb+buttsep,buttw,butth],...
          'string','Print',...
          'call','dfield6(''print'')','visible','off'); 
      
      dud.butt = dbutt;
      hhsetup = get(0,'showhiddenhandles');
      set(0,'showhiddenhandles','on');
 
     
      % File menu
      
      fmenu = findobj(gcf,'label','&File'); 
      delete(findobj(fmenu,'label','&New Figure'));
      delete(findobj(fmenu,'label','&Open...'));
      delete(findobj(fmenu,'label','&Close'));
      set(findobj(fmenu,'label','&Save'),...
          'separator','off',...
          'pos',1);
      set(findobj(fmenu,'label','Save &As...'),...
          'pos',2);
      delete(findobj(fmenu,'label','Pre&ferences...'));
      set(findobj(fmenu,'label','&Export...'),...
          'pos',3);
      set(findobj(fmenu,'label','Pa&ge Setup...'),'pos',4,...
                        'separator','on');
      set(findobj(fmenu,'label','Print Pre&view...'),'pos',5);
      set(findobj(fmenu,'label','Print Set&up...'),'pos',6);
      set(findobj(fmenu,'label','&Print...'),'pos',7);
      
      mequit = uimenu(fmenu,'label','Quit DFIELD6',...
          'call','dfield6(''quit'')','separator','on');
      
      % Tools menu
      tmenu = findobj(gcf,'label','&Tools'); 
      delete(tmenu);
            
% $$$       set(findobj(tmenu,'label','&Rotate 3D'),'vis','off');
% $$$       set(findobj(tmenu,'label','Zoom &Out'),'vis','off');
% $$$       set(findobj(tmenu,'label','&Zoom In'),'vis','off');
% $$$       set(findobj(tmenu,'label','Show Le&gend'),'vis','off');
% $$$       set(findobj(tmenu,'label','Enable &Plot Editing'),'vis','off');
% $$$       set(findobj(tmenu,'label','Show &Toolbar'),...
% $$$     'call','dfield6(''showbar'')')

      % Edit menu
      
      emenu = findobj(gcf,'label','&Edit'); 
      menu(2) = uimenu(emenu,'label','Zoom in.',...
                       'call','dfield6(''zoomin'')',...
                       'pos',1);
      
      zbmenu = uimenu(emenu,'label','Zoom back.',...
                      'call','dfield6(''zoomback'')',...
                      'enable','off',...
                      'tag','zbmenu',...
                      'pos',2);
      
      metext = uimenu(emenu,'label','Enter text on the Display Window.',...
                      'call','dfield6(''text'')',...
                      'separator','on',...
                      'pos',3);
         
      medall = uimenu(emenu,'label','Erase all solutions.',...
                      'call','dfield6(''dallsol'')',...
                      'separator','on',...
                      'pos',4);
      
      medalllev = uimenu(emenu,'label','Erase all level curves.',...
                         'call','dfield6(''dalllev'')',...
                         'separator','off');
      
      medall = uimenu(emenu,'label','Erase all graphics objects.',...
                      'call','dfield6(''dall'')',...
                      'separator','off',...
                      'pos',5);
         
      medel = uimenu(emenu,'label','Delete a graphics object.',...
                     'call','dfield6(''delete'')',...
                     'visible','on',...
                     'pos',6);
      set(findobj(emenu,'label','&Undo'),'separator','on',...
                        'pos',7);
      set(findobj(emenu,'label','Cu&t'),'pos',8);
      set(findobj(emenu,'label','&Copy'),'pos',9);
      set(findobj(emenu,'label','&Paste'),'pos',10);
      set(findobj(emenu,'label','&Select All'),'pos',11);
      set(findobj(emenu,'label','Copy &Figure'),'pos',12);
      set(findobj(emenu,'label','Copy &Options'),'pos',13);
      set(findobj(emenu,'label','F&igure Properties'),'pos',14);
      set(findobj(emenu,'label','&Axes Properties'),'pos',15);
      set(findobj(emenu,'label','C&urrent Object Properties'),'pos',16);

      % Options menu
            
      menu(1) = uimenu('label','Options','visible','off','pos',3);
      
      menukey = uimenu(menu(1),'label','Keyboard input.','call',...
          'dfield6(''kbd'')');
      
      mesev   = uimenu(menu(1),'label','Plot several solutions.',...
          'call','dfield6(''several'')');
      
      menulevel = uimenu(menu(1),'label','Plot level curves.',...
                         'call','dfield6(''level'')',...
                         'separator','off','tag','level');

      plcall = [
          'ud = get(gcf,''user'');',...
          'me = gcbo;',...
          'label = get(me,''label'');',...  
          'if strcmp(label,''Show the phase line.''),',...
          '  set(me,''label'',''Hide the phase line.'');',...
          '  ud.pline = ''on'';',...
          '  set(ud.plineh,''vis'',''on'');',...          
          '  set(gcf,''user'',ud);',...
          'else,',...
          '  set(me,''label'',''Show the phase line.'');',...
          '  ud.pline = ''off'';',...
          '  set(ud.plineh,''vis'',''off'');',...         
          '  kk = find(ishandle(ud.plhand));',...         
          '  delete(ud.plhand(kk));',...
          '  ud.plhand = [];',...         
          '  set(gcf,''user'',ud);',...
          'end'];  
      
      mepline = uimenu(menu(1),'label','Show the phase line.',...
          'call',plcall,'tag','pline','separator','on');          
      
      mesolvset = uimenu(menu(1),'label','Solver settings.',...
          'call','dfield6(''solvset'')');
      mesolve = uimenu(menu(1),'label','Solver.');
      
      solverstr = ['ud = get(gcf,''user'');',...
            'me = gcbo;',...
            'meud = get(me,''user'');',...
            'ud.settings.refine = meud.refine;',...
            'ud.settings.tol = meud.tol;',...
            'ud.settings.solver = meud.solver;',...
            'ud.settings.stepsize = meud.stepsize;',...
            'ud.settings.hmax = meud.hmax;',...
            'set(ud.solver,''checked'',''off'');',...
            'set(me,''checked'',''on'');',...
            'set(gcf,''user'',ud);',...
            'dfset = findobj(''name'',''DFIELD6 Setup'');',...
            'sud = get(dfset,''user'');',...
            'sud.settings = ud.settings;',...
            'set(dfset,''user'',sud);',...
            'dfield6(''solvset'');'];
      solver = dud.settings.solver;
      dpset.refine = 4;
      dpset.tol = 5e-4;
      dpset.solver = 'Dormand-Prince';
      dpset.stepsize = 0;
      dpset.hmax = 0;
      if strcmp(solver,'Dormand-Prince')
        dpch = 'on';
      else
        dpch = 'off';
      end
      
      eulset.refine = 1;
      eulset.tol = 0;
      eulset.solver = 'Euler';
      eulset.stepsize = 0.1;
      eulset.hmax = 0;
      if strcmp(solver,'Euler')
        eulch = 'on';
      else
        eulch = 'off';
      end
      rk2set.refine = 1;
      rk2set.tol = 0;
      rk2set.solver = 'Runge-Kutta 2';
      rk2set.stepsize = 0.1;
      rk2set.hmax = 0;
      if strcmp(solver,'Runge-Kutta 2')
        rk2ch = 'on';
      else
        rk2ch = 'off';
      end
      rk4set.refine = 1;
      rk4set.tol = 0;
      rk4set.solver = 'Runge-Kutta 4';
      rk4set.stepsize = 0.1;
      rk4set.hmax = 0;
      if strcmp(solver,'Runge-Kutta 4')
        rk4ch = 'on';
      else
        rk4ch = 'off';
      end
      
      dud.solver(1) = uimenu(mesolve,'label','Dormand-Prince',...
          'checked',dpch,...
          'call',solverstr,'user',dpset);
      
      dud.solver(2) = uimenu(mesolve,'label','Euler',...
          'checked',eulch,...
          'call',solverstr,'user',eulset);
      
      dud.solver(3) = uimenu(mesolve,'label','Runge-Kutta 2',...
          'checked',rk2ch,...
          'call',solverstr,'user',rk2set);
      
      dud.solver(4) = uimenu(mesolve,'label','Runge-Kutta 4',...
          'checked',rk4ch,...
          'call',solverstr,'user',rk4set); 
      
      medir = uimenu(menu(1),'label','Solution direction.');
      
      directionstr = ['ud = get(gcf,''user'');',...
            'me = gcbo;',...
            'ud.dir = get(me,''user'');',...
            'set(ud.direction,''checked'',''off'');',...
            'set(me,''checked'',''on'');',...
            'set(gcf,''user'',ud);'];
      
      dud.direction(1) = uimenu(medir,'label','Both',...
          'checked','on',...
          'user',0,...
          'call',directionstr);
      dud.dir = 0;
      
      dud.direction(2) = uimenu(medir,'label','Forward',...
          'user',1,...
          'call',directionstr);
      
      dud.direction(3) = uimenu(medir,'label','Back',...
          'user',-1,...
          'call',directionstr);
      
      meset   = uimenu(menu(1),'label','Window settings.',...
          'call','dfield6(''settings'')');
      
      menu(6) = uimenu(menu(1),'label','Make the Display Window inactive.',...
          'call','dfield6(''hotcold'')','separator','on');
      
      % View menu
      
      set(findobj(gcf,'label','&View'),'pos',4);
      set(findobj(gcf,'label','&Figure Toolbar'),...
          'call','dfield6(''showbar'')');
      dud.menu = menu;
      
      set(dfdisp,'ToolBar','none');
      set(0,'showhiddenhandles',hhsetup);
      set(gcf,'WindowButtonDownFcn','dfield6(''down'')');
      dud.uicont = [nframe,dud.notice,dbutt([1 2 3])];
      hh1 = [dud.axes,nframe,dud.notice,dbutt([1 2 3])];
      set(hh1,'units','norm');
      hh2 = [nframe,dud.notice,dbutt(2:3),dud.axes,dud.menu(1)];
      set(hh2,'visible','on');
      
      set(dfdisp,'vis','on');
      dud.printstr = 'print -noui';
      
      dud.pline = 'off';    
      dud.arr = [];
      dud.solhand = [];
      dud.plhand = [];
      dud.level = ' ';
      dud.contours = zeros(0,1);
    end   % if (~isempty(dfdisp))  & else
    
    dfdispa = dud.axes;
    axes(dfdispa);
    xlabel(Tname);
    ylabel(Xname);

    hv = get(0,'showhiddenhandles');
    set(0,'showhiddenhandles','on');
    set(findobj('tag','zbmenu'),'enable','off');
    set(0,'showhiddenhandles',hv);
    k = find(pflag);
    if ~isempty(k)
      lk = length(k);
      switch lk
        case 1
          tstr = [tstr,'             ',pstring{k}];
        case 2
          tstr = [tstr,'             ',pstring{1},'    ',pstring{2}];
      end
    end
    
    title(tstr)
    
    % Initialize important information as user data.
    
    dud.function = dfcn;
    if ~isempty(dud.solhand)    
      delete(dud.solhand);
    end    
    dud.solhand = [];   % Handles to solution curves.
    if isstruct(dud.arr)
      hand = [dud.arr.lines;dud.arr.arrows];
      delete(hand);
    end
    dud.arr = [];       % Handles for the direction and vector fields.
    dud.wmat = [];
    set(findobj('tag','pline'),'label','Show the phase line.');    
    dud.pline = 'off';
    kk = find(ishandle(dud.plhand));    
    delete(dud.plhand(kk));
    dud.plhand = [];    
   
    dud.color = sud.color;
    dud.fieldtype = Arrflag;
    dud.npts = NumbFPts;  

    ud.y = zeros(2,1);
    ud.i = 0;
    ud.line = 0;
    wind = dud.syst.wind(:);
    dwind = [wind(1); wind(3); -wind(2); -wind(4)];
    DY = [wind(2)-wind(1); wind(4)-wind(3)];
    ud.DY = DY;  
    ud.cwind = wind - dud.settings.magn*[DY(1);-DY(1);DY(2);-DY(2)];
    ud.stop = 0;
    ud.gstop = 1;
    ud.plot = 1;
    hmax = DY(1)/10;
    dud.settings.hmax = hmax;
    mud = get(dud.solver(1),'user');
    mud.hmax = hmax;
    set(dud.solver(1),'user',mud);
    set(dfdisp,'user',dud);
    set(dud.axes,'user',ud);
    ppkbd = findobj('name','DFIELD6 Keyboard input','vis','on');
    if ~isempty(ppkbd),dfield6('kbd'),end
    dfield6('dirfield',dfdisp); 
end

 
elseif strcmp(action,'dirfield')
  
  % 'dirfield' computes and plots the field elements.  This is the entry
  % point both from 'proceed' and from later commands that require the
  % recomputation of the field elements.
  
  % Find DFIELD6 Display and get the user data.
  
  disph = input1;   
  
  dud = get(disph,'user');
  color = dud.color;
  
  dfcn = dud.function;
  dfdispa = dud.axes;
  WINvect = dud.syst.wind;
  settings = dud.settings;
  notice = dud.notice;
  if notice
    nstr = get(notice,'string');
    nstr(1:4)=nstr(2:5);
    nstr{5,1} = 'Computing the field elements.';
    set(notice,'string',nstr);
    
    % Augment the window matrix
    
    wmat = dud.wmat;
    wrows = size(wmat,1);
    wflag = 0;
    for k = 1:wrows
      if wmat(k,:)==WINvect
        wflag = 1;
      end
    end
    if wflag == 0
      wmat = [wmat;WINvect];
      dud.wmat = wmat;
    end
    if wrows 
      hv = get(0,'showhiddenhandles');
      set(0,'showhiddenhandles','on');
      hhh = findobj('tag','zbmenu');
      set(hhh,'enable','on');
      set(0,'showhiddenhandles',hv);
    end
  end
  
  
  Tmin = WINvect(1);
  Tmax = WINvect(2);
  Xmin = WINvect(3);
  Xmax = WINvect(4);
  
  N = dud.npts;
  deltax=(Xmax - Xmin)/(N-1);
  deltat=(Tmax - Tmin)/(N-1);
  
  % Set up the display window.
  
  Dxint=[Xmin-deltax,Xmax+deltax];
  Dtint=[Tmin-deltat,Tmax+deltat];
  
  % Set up the original mesh.
  
  XXXg=Xmin + deltax*[0:N-1];
  TTTg=Tmin + deltat*[0:N-1];
  
  [Tt,Xx]=meshgrid(TTTg,XXXg);
  
  % Calculate the line and vector fields.
  
  Xx=Xx(:);Tt=Tt(:);
  Ww = zeros(size(Xx));
  Ww = feval(dfcn,Tt',Xx');
  Vv = ones(size(Ww)) + Ww*sqrt(-1);
  Vv = Vv.';
  Arrflag = dud.fieldtype;
  
  
  mgrid = Tt+Xx.*sqrt(-1); % mgrid = mgrid(:);
  zz=Vv.';
  sc = min(deltat,deltax);
  
  arrow=[-1,1].';
  zzz=sign(zz);
  scale = sqrt((real(zzz)/deltat).^2+(imag(zzz)/deltax).^2);
  ww = (zzz == 0);
  scale = scale + ww;
  aa1 = 0.3*arrow*(zzz./scale)+ones(size(arrow))*(mgrid.');
  [r,c] = size(aa1);
  aa1=[aa1;NaN*ones(1,c)];
  aa1=aa1(:);
  
  arrow = [0,1,.7,1,.7].' + [0,0,.25,0,-.25].' * sqrt(-1);
  zz=sign(zz).*((abs(zz)).^(1/3));
  scale = 0.9*sc./max(max(abs(zz)));
  aa2 = scale*arrow*zz +ones(size(arrow))*(mgrid.');
  [r,c] = size(aa2);
  aa2=[aa2;NaN*ones(1,c)];
  aa2=aa2(:);
  axes(dfdispa);
  
  arr = dud.arr;        % Delete the old field data.
  if isstruct(arr)
    hand = [arr.lines;arr.arrows];  
    delete(hand);
  end
  arrh1 = plot(real(aa1),imag(aa1),'color',color.arrows,'visible','off');
  arrh2 = plot(real(aa2),imag(aa2),'color',color.arrows,'visible','off');

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

  switch Arrflag
    case 'lines'
      set(arrh1,'visible','on');
    case 'arrows'
      set(arrh2,'visible','on');
  end
  dud.arr.lines = arrh1;           % Save the handles for later use.
  dud.arr.arrows = arrh2;          % Save the handles for later use.

  if notice
    nstr = get(notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5,1} = 'Ready.';
    set(notice,'string',nstr);
  end
  axis([Dtint,Dxint]);
  aa = Dtint(1) + (Dtint(2) - Dtint(1))/100;
  if ~isempty(dud.plhand)
    set(dud.plhand,'xdata',aa);
  end
  if (isfield(dud,'plineh') & (~isempty(dud.plineh)))
    set(dud.plineh,'xdata',[aa,aa],'ydata',Dxint);
  else
    dud.plineh = plot([aa,aa],Dxint,'color',dud.color.pline);
  end
  set([dud.plhand;dud.plineh],'vis',dud.pline);  
  set(disph,'user',dud);

elseif strcmp(action,'hotcold')

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

  dfdisp = gcf;
  dud = get(dfdisp,'user');
  nstr = get(dud.notice,'string');
  nstr(1:4) = nstr(2:5);
  mehc = dud.menu(6);
  if (findstr(get(mehc,'label'),'inactive'))
    set(dfdisp,'WindowButtonDownFcn',' ');
    set(mehc,'label','Make the Display Window active.');
    nstr{5,1} = 'The Display Window is not active.';
    set(dud.notice,'string',nstr);
  else
    set(dfdisp,'WindowButtonDownFcn','dfield6(''down'')');
    set(mehc,'label','Make the Display Window inactive.');
    nstr{5,1} = 'The Display Window is active.';
    set(dud.notice,'string',nstr);
  end
        
elseif strcmp(action,'down')
  
  % 'down' is the Window Button Down call.  It starts the computation of
  % solutions from a click of the mouse.
  
  disph = gcf;
  seltype = get(disph,'selectiontype');
  if strcmp(seltype,'alt')
    dfield6('zoom');
    return
  end
  dud = get(disph,'user');
  ax = dud.axes;
  ch = findobj('type','uicontrol','enable','on');
%  mh = findobj('type','uimenu','enable','on');
  set(ch,'enable','inactive');
%  set(mh,'enable','off');
  wbdf = get(disph,'WindowbuttonDownFcn');
  set(disph,'WindowbuttonDownFcn','');
  axes(ax);
  initpt = get(ax,'currentpoint');
  initpt = initpt(1,[1,2]);
  dfield6('solution',initpt,disph);
  set(disph,'WindowbuttonDownFcn',wbdf);
%  set([ch;mh],'enable','on');
  set(ch,'enable','on');
  notice = dud.notice;
  if notice
    nstr = get(notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5,1} = 'Ready.';
    set(notice,'string',nstr)
  end
  
elseif strcmp(action,'several')
  
  % 'several' allows the user to pick several initial points at once.
  % This is not needed in X-windows, but it is on the Macintosh.
  
  disph = gcf;
  ch = findobj('type','uicontrol','enable','on');
%  mh = findobj('type','uimenu','enable','on');
  set(ch,'enable','inactive');
%  set(mh,'enable','off')
  wbdf = get(disph,'WindowbuttonDownFcn');
  set(disph,'WindowbuttonDownFcn',wbdf);
  dud = get(disph,'user');
  notice = dud.notice;
  if notice
    nstr = get(notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5,1} = ['Pick initial points with the mouse.  ',...
        'Enter "Return" when finished.'];
    set(notice,'string',nstr)
  end
  [X,Y]=ginput;
  NN = length(X);
  for k = 1:NN
    initpt = [X(k),Y(k)];
    dfield6('solution',initpt,disph);
  end
  if notice
    nstr = get(notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5,1} = 'Ready.';
    set(notice,'string',nstr);
  end
%  set([ch;mh],'enable','on');
  set(ch,'enable','on');
  set(disph,'WindowbuttonDownFcn','dfield6(''down'')');
  

elseif strcmp(action,'solution')

  % 'solution' effects the computation and (erasemode == none) plotting of
  % solutions.  It also stores the data as appropriate.
  
  disph = input2;
  dud = get(disph,'user');
  tcol = dud.color.temp;
  pcol = dud.color.orb;
  notice = dud.notice;
  initpt = input1(:);
  dfcn = dud.function;
  dfdispa = dud.axes;
  aud=get(dfdispa,'user');  
  wind = aud.cwind;
  AA = wind(1);
  BB = wind(2);
  settings = dud.settings;
  ptstr = [' (',num2str(initpt(1),2), ', ', num2str(initpt(2),2), ')'];
  refine = settings.refine;
  tol = settings.tol;
  ud = get(dud.axes,'user');
  rtol = tol;
  atol = tol*ud.DY*1e-4';
  if length(initpt)  == 2
    switch dud.dir
     case 0
      intplus = [initpt(1), BB];
      intminus = [initpt(1), AA]; 
     case -1
      intplus = [initpt(1), initpt(1)];
      intminus = [initpt(1), AA]; 
     case 1
      intplus = [initpt(1), BB];
      intminus = [initpt(1), initpt(1)];
    end
  else
    intplus = [initpt(1),initpt(4)];
    intminus = [initpt(1),initpt(3)];
    initpt = initpt([1:2]);
  end       
  stopbutt =findobj('tag','stop');
  set(stopbutt,'vis','on','enable','on');

  cflag = 0;
  if intplus(2)>intplus(1)
    cflag = cflag + 1;
    if notice
      nstr = get(notice,'string');
      nstr(1:4) = nstr(2:5);
      nstr{5} = ['The forward orbit from',ptstr];
      set(notice,'string',nstr);
    end
    drawnow                             
    
    switch settings.solver
     case 'Dormand-Prince'
      [tp,xp] = dfdp45(dfcn,intplus,initpt(2),disph);
      
     case 'Euler'       
      [tp,xp] = dfeul(dfcn,intplus,initpt(2),disph);
      
     case 'Runge-Kutta 2'       
      [tp,xp] = dfrk2(dfcn,intplus,initpt(2),disph);
      
     case 'Runge-Kutta 4'       
      [tp,xp] = dfrk4(dfcn,intplus,initpt(2),disph);
      
     otherwise
      error('Undefined solver.');
    end
    aud = get(dfdispa,'user');
    hnew1 = aud.line;
    set(aud.pline,'erase','normal','color', pcol);  
    dud.plhand=[dud.plhand;aud.pline];   
  end           
  
  if intminus(2) < intminus(1)
    cflag = cflag + 2;
    if notice   
      nstr = get(notice,'string');
      nstr(1:4) = nstr(2:5);
      nstr{5} = ['The backward orbit from',ptstr];
      set(notice,'string',nstr);
    end
    drawnow
    
    switch settings.solver
     case 'Dormand-Prince'
      [tm,xm] = dfdp45(dfcn,intminus,initpt(2),disph);
      
     case 'Euler'       
      [tm,xm] = dfeul(dfcn,intminus,initpt(2),disph);
      
     case 'Runge-Kutta 2'       
      [tm,xm] = dfrk2(dfcn,intminus,initpt(2),disph);
      
     case 'Runge-Kutta 4'       
      [tm,xm] = dfrk4(dfcn,intminus,initpt(2),disph);
      
     otherwise
      error('Undefined solver.');
    end
    aud = get(dfdispa,'user');
    hnew2 = aud.line;
    set(aud.pline,'erase','normal','color', pcol);     
    dud.plhand=[dud.plhand;aud.pline];   
    set(stopbutt,'vis','off');
  end     % if intminus(2) < intminus(1)        
  
  % Store the trajectory.
  
  switch cflag
   case 1 % positive only
    set(hnew1,'xdata',tp(:),'ydata',xp(:),'color',pcol);
    set(hnew1,'erase','normal');
    dud.solhand = [dud.solhand;hnew1];
   case 2 % negative only
    x = flipud(xm);
    t = flipud(tm);
    set(hnew2,'xdata',t(:),'ydata',x(:),'color',pcol);
    set(hnew2,'erase','normal');
    dud.solhand = [dud.solhand;hnew2];
    
   case 3 % both directions
    x = flipud(xm);
    t = flipud(tm);
    x=[x;xp];
    t=[t;tp];
    delete(hnew1);
    set(hnew2,'xdata',t(:),'ydata',x(:),'color',pcol);
    set(hnew2,'erase','normal');
    dud.solhand = [dud.solhand;hnew2];
  end    % switch cflag
  set(disph,'user',dud);
  
elseif strcmp(action,'kcompute')
  
  % 'kcompute' is the call back for the Compute 
  % button on the DFIELD6 Keyboard figure.
  
  compute = 1;
  kh = get(gcf,'user');
  dfdisp = findobj('name','DFIELD6 Display');
  if (isempty(dfdisp))
    dfield6('confused');
  end
  dud = get(dfdisp,'user');
  dfdispa = dud.axes;
  aud = get(dfdispa,'user');
  ch = findobj('type','uicontrol','enable','on');
%  mh = findobj('type','uimenu','enable','on');
  set(ch,'enable','inactive');
%  set(mh,'enable','off')
  set(dfdisp,'WindowbuttonDownFcn','');
  x = str2num(get(kh.xval,'string'));
  tzero = str2num(get(kh.tval,'string'));
  if get(kh.spec,'value')
    t0 = str2num(get(kh.t0,'string'));
    tf = str2num(get(kh.tf,'string'));
    initpt = [tzero,x,t0,tf];
    if (length(initpt) ~= 4)
      warndlg({'Values must be entered for all of the entries.'},...
          'Illegal input');
      compute = 0;
    elseif tf <= t0
      warndlg({'The final time of the computation interval';...
      'must be smaller than the initial time.'},'Illegal input');
      compute = 0;
    elseif (tzero < t0) | (tzero > tf)
      warndlg('The initial time must be in the computation interval.',...
          'Illegal input');
      compute = 0;
    end
    aud.gstop = 0;
    set(dfdispa,'user',aud);
  else
    initpt = [tzero,x];
    if (length(initpt) ~= 2)
      warndlg({'Values must be entered for all of the entries.'},...
          'Illegal input');
      compute = 0;
    end
  end   % if get(kh.spec,'value')
  
  if compute
    dfield6('solution',initpt,dfdisp);
  end
  nstr = get(dud.notice,'string');
  nstr(1:4) = nstr(2:5);
  nstr{5} = 'Ready.';
  set(dud.notice,'string',nstr)
%  set([ch;mh],'enable','on');
  set(ch,'enable','on');
  set(dfdisp,'WindowbuttonDownFcn','dfield6(''down'')');
  aud.gstop = 1;
  set(dfdispa,'user',aud);
  
  
elseif strcmp(action,'kbd')
  
  % 'kbd' is the callback for the Keyboard Input menu selection.  It 
  % sets up the DFIELD6 Keyboard figure which allows accurate input of 
  % initial values using the keyboard.
  
  ppkbd = findobj('name','DFIELD6 Keyboard input');
  if ~isempty(ppkbd)
    delete(ppkbd);
  end    
  dfdisp = gcf;
  dud = get(dfdisp,'user');
  Xname = dud.syst.xname;
  Tname = dud.syst.tname;
  xnstr = ['The initial value of ',Xname,' = '];
  tnstr = ['The initial value of ',Tname,' = '];
  
  ppkbd = figure('name','DFIELD6 Keyboard input',...
      'vis','off',...
      'numb','off','tag','DFIELD6');
  
  dfield6('figdefault',ppkbd);
  set(ppkbd,'menubar','none');
  
  kbd.fr1 = uicontrol('style','frame');
  
  kbd.fr2 = uicontrol('style','frame');
  
  kbd.inst = uicontrol('style','text','horiz','center',...
      'string','Enter the initial conditions:');
  
  kbd.xname = uicontrol('style','text',...
      'horiz','right','string',xnstr);
  
  kbd.tname = uicontrol('style','text',...
      'horiz','right',...
      'string',tnstr);
  
  kbd.tval = uicontrol('style','edit');
  
  kbd.xval = uicontrol('style','edit',...
      'string','','call','');
  
  kbd.t0 = uicontrol('style','edit');
  
  kbd.tf = uicontrol('style','edit');
  
  tstring = ['<= ',Tname,' <='];
  kbd.t = uicontrol('style','text','string',tstring);
  
  kbd.spec = uicontrol('style','check','horiz','center',...
      'string','Specify a computation interval.');
  
  kbd.comp = uicontrol('style','push',...               
      'string','Compute','call','dfield6(''kcompute'')');
  
  kbd.close = uicontrol('style','push',...
      'string','Close','call','close');
  
  nudge = 3; 
  left = 3; varl = 80; buttw = 60;
  xex = get(kbd.xname,'extent');
  ht = xex(4)+nudge;
  yex = get(kbd.tname,'extent');
  nwdth = max(xex(3),yex(3)) + nudge;
  varl = varl*ht/19;
  fr1bot = left + ht;
  fr1ht = 3*nudge + 2*ht;
  fr2bot = fr1bot + fr1ht;
  frw = 2*nudge + nwdth + varl;
  intb = fr1bot + nudge;
  specb = intb + ht + nudge;
  intext = get(kbd.t,'extent');
  tw = intext(3);
  newwidth = tw + 2*varl +2*nudge;
  frw = max(frw,newwidth);
  margin = (frw -tw - 2*varl)/2;
  t0l = left+margin;
  tl = t0l + varl;
  tfl = tl + tw;
  xnb = fr2bot + nudge;
  tnb = xnb + ht;
  nl = left+nudge;
  xnwind = [nl,xnb-.15*ht,nwdth,ht];
  tnwind = [nl,tnb-.15*ht,nwdth,ht];
  vl = nl + nwdth;
  xvwind = [vl,xnb,varl,ht];
  tvwind = [vl,tnb,varl,ht];
  instb = tnb + ht + nudge;
  instw = nwdth + varl;
  instwind = [nl,instb,instw,ht];
  fr2ht = 3*nudge + 3*ht;
  % frw = 2*nudge + nwdth + varl;
  fr2wind = [left,fr2bot,frw,fr2ht];
  fr1wind = [left,fr1bot,frw,fr1ht];
  specw = frw - 2*nudge;
  specwind = [nl,specb,specw,ht];
  t0wind = [t0l,intb,varl,ht];
  twind = [tl,intb-.15*ht,tw,ht];
  tfwind = [tfl,intb,varl,ht];
  figw = 2*left + frw;
  fight = 3*left + ht + fr1ht + fr2ht;
  figwind = [30,300,figw,fight];
  buttw = frw/2;
  sep = (figw - 2*buttw)/3;
  closel = left;
  compl = closel+buttw;
  clwind = [closel,left,buttw,ht];
  compwind = [compl,left,buttw,ht];
  set(ppkbd,'pos',figwind);
  set(kbd.fr1,'pos',fr1wind);
  set(kbd.fr2,'pos',fr2wind);
  set(kbd.inst,'pos',instwind);
  set(kbd.xname,'pos',xnwind);
  set(kbd.tname,'pos',tnwind);
  set(kbd.xval,'pos',xvwind);
  set(kbd.tval,'pos',tvwind);
  set(kbd.comp,'pos',compwind);
  set(kbd.close,'pos',clwind);
  set(kbd.tname,'pos',tnwind);
  set(kbd.tval,'pos',tvwind);
  set(kbd.spec,'pos',specwind);
  set(kbd.t0,'pos',t0wind);
  set(kbd.t,'pos',twind);
  set(kbd.tf,'pos',tfwind);
  speccall = [
    'ud = get(gcf,''user'');',...
        'if get(gcbo,''value''),',...
        '  set([ud.t0,ud.t,ud.tf],''vis'',''on'');',...
        'else,',...
        '  set([ud.t0,ud.t,ud.tf],''vis'',''off'');',...
        'end'];
  
  set(kbd.spec,'call',speccall);
  set(ppkbd,'resize','on');
  set(findobj(ppkbd,'type','uicontrol'),'units','normal');
  
  set(kbd.spec,'value',0);
  set(ppkbd,'user',kbd,'vis','on');
  set([kbd.t0,kbd.t,kbd.tf],'vis','off')
  set(findobj(ppkbd,'type','uicontrol'),'units','normal');
  figure(ppkbd)
  
elseif strcmp(action,'zoomin')
  
  % 'zoomin' is the callback for the Zoomin menu item.  It allows the
  % user to pick a new display rectangle.
  
  set(gcf,'WindowButtonDownFcn','dfield6(''zoom'')',...
      'WindowButtonUpFcn','1;','inter','on');
  set(gca,'inter','on');
  dud = get(gcf,'user');
  nstr = get(dud.notice,'string');
  nstr(1:4) = nstr(2:5);
  nstr{5} = ['Pick a new display rectangle by clicking and ',...
      'dragging the mouse, or by clicking on a point.'];
  set(dud.notice,'string',nstr);
  
elseif strcmp(action,'zoom')
  
  disph = gcf;
  dud = get(disph,'user');
  axh = dud.axes;
  aud = get(axh,'user');
  DY = aud.DY;
  w = dud.syst.wind;
  q1 = get(disph,'currentpoint');
  p1 = get(axh,'currentpoint');
  p1 = p1(1,1:2);
  rbbox([q1 0 0],q1);
  p2 = get(axh,'currentpoint');
  p2 = p2(1,1:2);
  if all(abs(p1'-p2')>0.01*DY)
    a = [p1;p2];
    a = [min(a);max(a)];
    DY = (a(2,:) - a(1,:))';  
  else
    DY = DY/4;
    a(1) = max(w(1),p1(1)-DY(1));
    a(2) = min(w(2),p1(1)+DY(1));
    a(3) = max(w(3),p1(2)-DY(2));
    a(4) = min(w(4),p1(2)+DY(2));
    DY(1) = a(2) - a(1);
    DY(2) = a(4) - a(3);
  end    
  WINvect = a(:)';
  dud.syst.wind = WINvect;
  aud.DY = DY;
  dwind = [WINvect(1); WINvect(3); -WINvect(2); -WINvect(4)];
  aud.cwind = WINvect(:) - ...
      dud.settings.magn*[DY(1);-DY(1);DY(2);-DY(2)];
  set(axh,'user',aud);  
  set(disph,'user',dud);
  set(disph,'WindowButtonDownFcn','dfield6(''down'')',...
      'WindowButtonUpFcn','');
  dfield6('dirfield',disph);
  dfset = findobj('name','DFIELD6 Setup');
  if isempty(dfset)
    dfield6('confused');
  else
    sud = get(dfset,'user');
    sud.c.wind = WINvect;
    sud.o.wind = WINvect;
    set(sud.h.wind(1),'string',num2str(WINvect(1)));
    set(sud.h.wind(2),'string',num2str(WINvect(2)));
    set(sud.h.wind(3),'string',num2str(WINvect(3)));
    set(sud.h.wind(4),'string',num2str(WINvect(4)));
    set(dfset,'user',sud);
  end
  
elseif strcmp(action,'dall')
  
  % 'dall' is the callback for the Erase all graphics objects.  
   
  disph = gcf;
  dud = get(disph,'user');
  notice = dud.notice;
  kk = find(ishandle(dud.solhand));    
  delete(dud.solhand(kk));
  dud.solhand = [];
  kk = find(ishandle(dud.plhand));    
  delete(dud.plhand(kk));
  dud.plhand = [];
  kk = find(ishandle(dud.contours));    
  delete(dud.contours(kk));
  dud.contours = [];
  if notice
    set(dud.butt(1),'enable','off');
  end
  set(disph,'user',dud);
  
elseif strcmp(action,'dallsol')
   
  % 'dallsol' is the callback for the Erase all solutions option.  
  
  disph = gcf;
  dud = get(disph,'user');
  notice = dud.notice;
  kk = find(ishandle(dud.solhand));      
  delete(dud.solhand(kk));
  dud.solhand = [];
  kk = find(ishandle(dud.plhand));    
  delete(dud.plhand(kk));
  dud.plhand = [];
  if notice
    set(dud.butt(1),'enable','off');
  end
  set(disph,'user',dud);
  
elseif strcmp(action,'dalllev')
   
  % 'dalllev' is the callback for the Erase all level curves option.  
  
  disph = gcf;
  dud = get(disph,'user');
  notice = dud.notice;
  kk = find(ishandle(dud.contours));       
  delete(dud.contours(kk));
  dud.contours = [];
  if notice
    set(dud.butt(1),'enable','off');
  end
  set(disph,'user',dud);
   
elseif strcmp(action,'solvset')

  dud = get(gcf,'user');
  name = dud.settings.solver;
  
  setfig = findobj('name','DFIELD6 Solver settings');
  if ~isempty(setfig)
    data = get(setfig,'user');
    if strcmp(data.settings.solver,name)
      figure(setfig);
      return
    else
      close(setfig);
    end
  end
  data.settings = dud.settings;
      
  nstepcall =['data = get(gcf,''user'');',...
        'ss = max(round(str2num(get(data.d1,''string''))),1);',...
        'data.settings.refine = ss;',...
        'set(data.d1,''string'',num2str(ss));',...
        'set(data.chb,''enable'',''on''),',...
        'set(gcf,''user'',data);'];

  ssizecall = ['data = get(gcf,''user'');',...
        'data.settings.stepsize = str2num(get(data.d1,''string''));',...
        'set(data.d1,''string'',num2str(data.settings.stepsize));',...
        'set(data.chb,''enable'',''on''),',...
        'set(gcf,''user'',data);'];
  
  chcall = ['data = get(gcf,''user'');',...
        'dfdisp = findobj(''name'',''DFIELD6 Display'');',...
        'dud = get(dfdisp,''user'');',...
        'dud.settings = data.settings;',...
        'set(dfdisp,''user'',dud);',...
        'dfset = findobj(''name'',''DFIELD6 Setup'');',...
        'ud = get(dfset,''user'');',...
        'ud.settings = data.settings;',...
        'set(dfset,''user'',ud);',...
        'set(data.chb,''enable'',''off'');'];
  
  speedcall = ['data = get(gcf,''user'');',...
        'me = data.speed;',...
        'val = round(get(me,''value''));',...
        'set(me,''value'',val);',...
        'set(data.sp.val,''string'',num2str(val));',...
        'data.settings.speed = val;',...
        'set(data.chb,''enable'',''on''),',...
        'set(gcf,''user'',data);'];
  
  switch name
  case 'Dormand-Prince'
    ssname = 'Settings for the Dormand-Prince solver.';
    d1call = nstepcall;
    d1string = 'Number of plot steps per computation step:  ';
    d1data = num2str(data.settings.refine);
  case 'Euler'
    ssname = 'Settings for Euler''s method.';
    d1call = ssizecall;
    d1string = 'Step size:     ';
    d1data = num2str(data.settings.stepsize);
  case 'Runge-Kutta 2'
    ssname = 'Settings for the second order Runge-Kutta method.';
    d1call = ssizecall;
    d1string = 'Step size:     ';
    d1data = num2str(data.settings.stepsize);
  case 'Runge-Kutta 4'
    ssname = 'Settings for the fourth order Runge-Kutta method.';
    d1call = ssizecall;
    d1string = 'Step size:     ';
    d1data = num2str(data.settings.stepsize);
  end

  left = 5; nudge = 1; varl = 70;
        
  setfig = figure('name','DFIELD6 Solver settings','numb','off',...
      'tag','DFIELD6','vis','off');
  
  dfield6('figdefault',setfig);
  set(setfig,'menubar','none');

  setfr = uicontrol('style','frame');   
  
  setspfr = uicontrol('style','frame'); 
  
  d1 = uicontrol('style','text','horiz','left',...
      'string',d1string);
      
  data.d1 = uicontrol('style','edit','string',d1data,...
      'call',d1call);

  rtolcall =['data = get(gcf,''user'');',...
        'data.settings.tol = str2num(get(data.rtol,''string''));',...
        'set(data.rtol,''string'',num2str(data.settings.tol));',...
        'set(data.chb,''enable'',''on''),',...
        'set(gcf,''user'',data);'];
  
  hmaxcall = ['data = get(gcf,''user'');',...
        'data.settings.hmax = str2num(get(data.hmax,''string''));',...
        'set(data.hmax,''string'',num2str(data.settings.hmax));',...
        'set(data.chb,''enable'',''on''),',...
        'set(gcf,''user'',data);'];
 
  
  ss = uicontrol('style','text','horiz','center',...
      'string',ssname);

  gob = uicontrol('style','push','string','Go Away','call','close');
  
  data.chb = uicontrol('style','push',...
      'string','Change settings',...
      'call',chcall,...
      'enable','off');
  
  data.speed = uicontrol('style','slider',...
      'string','Steps per second.',...
      'min',2,...
      'max',100,...
      'value',data.settings.speed,...
      'sliderstep',[1/98,10/98],...
      'call',speedcall);
  
  data.sp.min = uicontrol('style','text','string',' 2',...
      'horiz','left');
  data.sp.max = uicontrol('style','text','string','100 ',...
      'horiz','right');
  data.sp.val = uicontrol('style','text',...
      'string',num2str(data.settings.speed),...
      'horiz','center');
  
  pps = uicontrol('style','text',...
      'string','Solution steps per second.');
  ext1 = get(d1,'extent');
  ht = ext1(4)+3;
  ext2 = get(ss,'extent');
  stw = ext1(3)+5;              
  varl = varl*ht/19;
  setspfrb = 2*left + 2*ht;
  setspfrht = 3*nudge + 3*ht;
  setfrb = setspfrb + setspfrht;
  setfrh = 4*nudge + 4*ht;
  figh = setfrb + setfrh + left;
  gobb = left;
  chbb = gobb + ht;
  hmaxb = setfrb + nudge;
  d1b = hmaxb + ht;
  rtolb = d1b + ht;
  ssizeb = (rtolb+d1b)/2;
  ssb = rtolb + ht + nudge;
  ssw = max(max(stw+varl,ext2(3)),300);
  frw = 2*nudge + ssw;
  figw = frw + 2*left;
  buttw = figw/3;
  sep = figw/9;
  sl = left + nudge;
  d1l = (figw -stw-varl)/2;
  gobl = sep;
  chbl = 2*sep + buttw;
  spb1 = setspfrb+2*nudge;
  spb2 = spb1+ht-nudge;
  spb3 = spb2+ht+nudge;
  sptw = varl;
  sptl1 = sl;
  sptsep = (ssw-3*sptw-sl)/2;
  sptl2 = sptl1 + sptw + sptsep;
  sptl3 = sptl2 + sptw + sptsep;

  sunit = get(0,'units');
  set(0,'units','pix');
  ssize = get(0,'screensize');
  figb = ssize(4) - figh - 50;
  
  set(setfig,'pos',[20,figb,figw,figh]);
  set(setspfr,'pos',[left,setspfrb,frw,setspfrht]);
  set(data.speed,'pos',[sl,spb1,ssw,ht],...
                 'backgroundcolor',0.6*[1 1 1],...
                 'foregroundcolor','r');
  set(data.sp.min,'pos',[sptl1,spb2,sptw,ht]);
  set(data.sp.max,'pos',[sptl3,spb2,sptw,ht]);
  set(data.sp.val,'pos',[sptl2,spb2,sptw,ht]);
  set(pps,'pos',[sl,spb3,ssw,ht]);
  set(setfr,'pos',[left,setfrb,frw,setfrh]);
  set(ss,'pos',[sl,ssb,ssw,ht]);
  set(d1,'pos',[d1l,d1b,stw,ht]);
  set(data.d1,'pos',[d1l+stw,d1b,varl,ht]);
  set(gob,'pos',[left,gobb,frw,ht]);
  set(data.chb,'pos',[left,chbb,frw,ht]);
  if strcmp(name,'Dormand-Prince')
    rtol = uicontrol('style','text','horiz','left',...
        'pos',[sl+5,rtolb,stw-5,ht],...
        'string','Relative error tolerance: '); 
    data.rtol = uicontrol('style','edit','call',rtolcall,...
        'pos',[sl+stw,rtolb,varl,ht],...
        'string',num2str(data.settings.tol));
    hmax = uicontrol('style','text','horiz','left',...
        'pos',[sl+5,hmaxb,stw-5,ht],...
        'string','Maximum step size: '); 
    data.hmax = uicontrol('style','edit','call',hmaxcall,...
        'pos',[sl+stw,hmaxb,varl,ht],...
        'string',num2str(data.settings.hmax));
    set(d1,'pos',[sl+5,d1b,stw-5,ht]);
    set(data.d1,'pos',[sl+stw,d1b,varl,ht]);
  end

  set(setfig,'user',data);
  set(setfig,'units','normal');
  set(findobj(setfig,'type','uicontrol'),'units','normal');
  set(setfig,'vis','on','resize','on');
  set(get(setfig,'child'),'vis','on');

elseif strcmp(action,'settings')
  
  % 'settings' is the call back for the Settings menu option.  It sets
  % up the DFIELD6 Settings window, which allows the user to interactively
  % change several parameters that govern the behaviour of the program.
        
  dud = get(gcf,'user');
  data.settings = dud.settings;
  fieldtype = dud.fieldtype;
  npts = dud.npts;
  data.fieldtype = fieldtype;
  data.npts = npts;
  data.magn = dud.settings.magn;
  setfig = findobj('name','DFIELD6 Window settings');
  if ~isempty(setfig)
    close(setfig);
  end
  setfig = figure('name','DFIELD6 Window settings',...
                  'numb','off',...
                  'tag','DFIELD6','vis','off');
  
  dfield6('figdefault',setfig);
  set(setfig,'menubar','none');

  dirffr = uicontrol('style','frame');  
  cwfr = uicontrol('style','frame');

  ss = uicontrol('style','text','horiz','center',...
      'string','The direction field.');
    
  kk = 1+2*data.settings.magn;
  magcall = ['data = get(gcf,''user'');',...
        'mag =  (str2num(get(data.mag,''string''))-1)/2;',...
        'data.magn = max(0,mag);',...
        'set(gcf,''user'',data);'];
  
  cw1 = uicontrol('style','text','horiz','left',...
      'string','The calculation window is');
  data.mag = uicontrol('style','edit','call',magcall,...
      'string',num2str(kk));
  cw2 = uicontrol('style','text','horiz','left',...
      'string',' times as large as the ');
  cw3 = uicontrol('style','text','horiz','left',...
      'string','display window.');
  
  gob = uicontrol('style','push','string','Go Away','call','close');
  
  chb = uicontrol('style','push',...
      'string','Change settings',...
      'call','dfield6(''setchange'')');

  left = 1; nudge = 3; 
  
  cw1ext = get(cw1,'extent');
  cw1w = cw1ext(3);
  cw2ext = get(cw2,'extent');
  cw2w = cw2ext(3);
  ht = cw1ext(4)+nudge;
 
  cwfrb = 2*left + ht;
  cwfrh = 4*nudge + 2*ht;
  dirffrb = cwfrb + cwfrh;
  rbfrh = 3*ht + 3*nudge;
  dirffrh = 4*nudge + 4*ht;
  figh = dirffrb + dirffrh + left;
  bb = left;
  cw3b = cwfrb + nudge; 
  cw2b = cw3b + ht;     % = cw1b = cweb
  rb1b = dirffrb + nudge;
  rb2b = rb1b + ht;
  rb3b = rb2b + ht;

  ssb = rb3b + ht + nudge;
  cwew = 40*ht/19;
  cww = cw1w + cwew + cw2w;

  frw = 2*nudge + cww;
  figw = frw + 2*left;
  buttw = frw/2;

  sl = left + nudge;
  gobl = left;
  chbl = left + buttw;
  
  sunit = get(0,'units');
  set(0,'units','pix');
  ssize = get(0,'screensize');
  figb = ssize(4) - figh - 50;
  
  set(setfig,'pos',[20,figb,figw,figh]);
  set(dirffr,'pos',[left,dirffrb,frw,dirffrh]);
  set(cwfr,'pos',[left,cwfrb,frw,cwfrh]);
  set(ss,'pos',[sl,ssb,cww,ht]);
  set(cw1,'pos',[sl,cw2b,cw1w,ht]);
  set(cw2,'pos',[sl+cw1w+cwew,cw2b,cw2w,ht]);
  set(cw3,'pos',[sl,cw3b,cww,ht]);
  set(data.mag,'pos',[sl+cw1w,cw2b,cwew,ht]);
  set(gob,'pos',[gobl,bb,buttw,ht]);
  set(chb,'pos',[chbl,bb,buttw,ht]);

  radw = .4*cww;
  rbl = sl + nudge;
  rbw = radw -4*nudge;
  typewindw = radw + 2*nudge;
  typewind = [left, dirffrb, typewindw, rbfrh];
  textwindl = left+typewindw;
  textleft = textwindl + nudge;
  textw = frw - typewindw;
  textwind = [textwindl, dirffrb,textw, rbfrh];
  typeframe = uicontrol('style','frame','pos',typewind,'visible','off');
  
  textframe = uicontrol('style','frame','pos',textwind,'visible','off');
  
  switch fieldtype
    case 'arrows'
      rval1 = 1;rval2 = 0;rval3 = 0;
    case 'lines'
      rval1 = 0;rval2 = 2;rval3 = 0;
    case 'none' 
      rval1 = 0;rval2 = 0;rval3 = 3;
    otherwise
      error(['Unknown fieldtype ',ud.o.fieldtype,'.'])
  end           

  data.rad(1) = uicontrol('style','radio',...
      'pos',[rbl rb3b rbw ht],...
      'string','Arrows','value',rval1,'visible','off');
  
  data.rad(2) = uicontrol('style','radio',...
      'pos',[rbl rb2b rbw ht],...
      'string','Lines','value',rval2,'max',2,'visible','off');
  
  data.rad(3) = uicontrol('style','radio',...
      'pos',[rbl rb1b rbw ht],...
      'string','None','value',rval3,'max',3,'visible','off');
 
  for i=1:3
    set(data.rad(i),'user',data.rad(:,[1:(i-1),(i+1):3]));
  end
  
  callrad = ['me = gcbo;',...
        'kk = get(me,''max'');',...
        'set(get(me,''user''),''value'',0),',...
        'set(me,''value'',kk);',...
        'ud = get(gcf,''user'');',...
        'switch kk,',...
        '       case 1, ud.fieldtype = ''arrows'';',...
        '       case 2, ud.fieldtype = ''lines'';',...
        '       case 3, ud.fieldtype = ''none'';',...
        'end,',...
        'set(gcf,''user'',ud);'];
  
  set(data.rad,'call',callrad);
 
  nfptsstr = {'Number of field points'; 'per row or column.'};
  nfptstext = uicontrol('style','text',...
      'pos',[textleft rb2b textw-2*nudge 1.5*ht],...
      'string',nfptsstr,'horizon','center','visible','off');
  
  callnfpts = ['me = gcbo;',...
        'kk = str2num(get(me,''string''));',...
        'if isempty(kk),',...
        '  set(me,''string'',''?'');',...
        '  kk = NaN;',...
        'else,',...
        '  kk = floor(kk);',...
        '  [m,N] = computer;'...
        '  if (N <= 8192),',...
        '       N = 32;',...
        '  else,',...
        '       N = 50;',...
        '  end,'...
        '  kk = min([N,max([5,kk])]);'...
        '  set(me,''string'',num2str(kk));'...
        'end,'...
  'data = get(gcf,''user'');'...
  'data.npts = kk;',...
  'set(gcf,''user'',data)'];
  
  npts = uicontrol('style','edit',...
      'pos',[textleft+(textw -30*ht/19)/2,rb1b 30*ht/19,ht],...
      'string',npts,'call',callnfpts,'visible','off');
   
  
  set(setfig,'user',data);
  set(setfig,'units','normal');
  set(findobj(setfig,'type','uicontrol'),'units','normal');
  set(setfig,'vis','on','resize','on');
  set(get(setfig,'child'),'vis','on');


elseif strcmp(action,'setchange')

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

  dfsett = gcf;
  data = get(dfsett,'user');
  dfdisp = findobj('name','DFIELD6 Display');
  dfset = findobj('name','DFIELD6 Setup');
  if isempty(dfdisp)|isempty(dfset)
    dfield6('confused');
    return
  end
  dud = get(dfdisp,'user');
  sud = get(dfset,'user');
  dud.settings.magn = data.magn;
  sud.settings.magn = data.magn;
  WINvect = dud.syst.wind;
  aud = get(dud.axes,'user'); 
  DY = aud.DY;    
  aud.cwind = WINvect(:) - dud.settings.magn*[DY(1);-DY(1);DY(2);-DY(2)];
  set(dfset,'user',sud);
  set(dfdisp,'user',dud);
  set(dud.axes,'user',aud);
  if dud.npts ~= data.npts
    dud.fieldtype = data.fieldtype;
    dud.npts = data.npts;
    sud.fieldtype = data.fieldtype;
    sud.npts = data.npts;
    set(dfdisp,'user',dud);
    set(dfset,'user',sud);
    dfield6('dirfield',dfdisp);
  elseif ~strcmp(dud.fieldtype,data.fieldtype)
    dud.fieldtype = data.fieldtype;
    sud.fieldtype = data.fieldtype;
    switch dud.fieldtype
    case 'lines'
      set(dud.arr.lines,'vis','on');
      set(dud.arr.arrows,'vis','off');
    case 'arrows'
      set(dud.arr.lines,'vis','off');
      set(dud.arr.arrows,'vis','on');
    case 'none'
      set(dud.arr.lines,'vis','off');
      set(dud.arr.arrows,'vis','off');
    end
    set(dfset,'user',sud);
    set(dfdisp,'user',dud);
  end
  close(dfsett)
  
elseif strcmp(action,'delete')

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

  disph = gcf;
  dud = get(disph,'user');
  arr = dud.arr;
  lv = get(arr.lines,'vis');
  av = get(arr.arrows,'vis');
  pv = get(dud.plineh,'vis');  
  handles = [arr.lines;arr.arrows;dud.plineh];
  set(handles,'vis','off');
  oldcall = get(disph,'WindowButtonDownFcn');
  set(disph,'WindowButtonDownFcn','');
  trjh = [dud.solhand;dud.plhand;dud.contours];
  notice = dud.notice;
  if notice
    nstr = get(notice,'string');
    nstr(1:4) = nstr(2:5);
    nstr{5} = 'Select a graphics object with the mouse.';
    set(notice,'string',nstr);
  end
  ginput(1);
  objh = get(disph,'currentobject');
  typ = get(objh,'type');
  axh = dud.axes;
  hh = [get(axh,'title'),get(axh,'xlabel'),get(axh,'ylabel')];
  if strcmp(typ,'line')
    dud.solhand = setdiff(dud.solhand,objh);  
    dud.plhand = setdiff(dud.plhand,objh);  
    dud.contours = setdiff(dud.contours,objh);  
    delete(objh); 
    if notice
      nstr(1:4) = nstr(2:5);
      nstr{5} = 'Ready.';
      set(notice,'string',nstr);
    end
  elseif strcmp(typ,'text') & ~ismember(objh,hh)
    delete(objh); 
    if notice
      nstr(1:4) = nstr(2:5);
      nstr{5} = 'Ready.';
      set(notice,'string',nstr);
    end
  else
    if notice
      nstr(1:4) = nstr(2:5);
      nstr{5} = 'The object you selected cannot be deleted.';
      set(notice,'string',nstr);
    end
  end

  set(arr.lines,'vis',lv);
  set(arr.arrows,'vis',av);
  set(dud.plineh,'vis',pv);  
  set(disph,'user',dud);
  set(disph,'WindowButtonDownFcn','dfield6(''down'')');
  
  
elseif strcmp(action,'print')
  
  dud = get(gcf,'user');
  nstr = get(dud.notice,'string');
  nstr(1:4) = nstr(2:5);
  nstr{5} = 'Preparing to print the DFIELD6 Display Window.  Please be patient.';
  set(dud.notice,'string',nstr);
  
  nstr(1:4) = nstr(2:5);
  nstr{5} = 'Printing the DFIELD6 Display Window.';
  set(dud.notice,'string',nstr);
  set(gcf,'pointer','watch');
  eval(dud.printstr);
  nstr(1:4) = nstr(2:5);
  nstr{5} = 'Ready.';
  set(dud.notice,'string',nstr);
  set(gcf,'pointer','arrow');
  
elseif strcmp (action,'zoomback')
  
  disph = gcf;
  dud = get(disph,'user');
  axh = dud.axes;
  Xname = dud.syst.xname;
  Tname = dud.syst.tname;
  wmat = dud.wmat;
  WINvect = dud.syst.wind;
  
  NN = size(wmat,1);
  
  wch = 0;j=0;
  while wch == 0
    j = j+1;
    if WINvect == wmat(j,:)
      wch = j;
    end
  end
  winstr = cell(1,NN);
  for j = 1:NN
    a = num2str(wmat(j,1));
    b = num2str(wmat(j,2));
    c = num2str(wmat(j,3));
    d = num2str(wmat(j,4));
    winstr{j} = [' ',a,' < ',Tname,' < ',b,'   &   ',c,' < ', Xname,' < ',d];
  end
  
  [sel,ok] = listdlg('liststring',winstr,...
     'selectionmode','single',...
     'listsize',[270,100],...
     'initialvalue',wch,...
     'name','DFIELD6 Zoomback',...
     'promptstring','Select a rectangle:',...
     'OKString','Zoom');
  
  if (~isempty(sel))
    WINvect = wmat(sel,:);
    dud.syst.wind = WINvect;
    set(gcf,'user',dud);
    dfset = findobj('name','DFIELD6 Setup');
    sud = get(dfset,'user');
    set(disph,'user',dud);
    set(sud.h.wind(1),'string',num2str(WINvect(1)));
    set(sud.h.wind(2),'string',num2str(WINvect(2)));
    set(sud.h.wind(3),'string',num2str(WINvect(3)));
    set(sud.h.wind(4),'string',num2str(WINvect(4)));
    sud.c.wind = WINvect;
    sud.o.wind = WINvect;
    set(dfset,'user',sud);
    aud = get(axh,'user');
    DY = [WINvect(2) - WINvect(1);WINvect(4) - WINvect(3)];
    aud.DY = DY;    
    aud.cwind = WINvect(:) - ...
             dud.settings.magn*[DY(1);-DY(1);DY(2);-DY(2)];
    set(axh,'user',aud);
    dfield6('dirfield',disph);
  end
  
elseif strcmp(action,'level')
  
  dfdisp = gcf;
  dud = get(dfdisp,'user');
  lfcn = dud.level;
  tname = dud.syst.tname;
  xname = dud.syst.xname;
  tname(find(abs(tname)==92))=[];  % Remove \s if any.
  xname(find(abs(xname)==92))=[];
  dflevel = findobj('name','DFIELD6 Level sets');
  if ~isempty(dflevel)
    delete(dflevel)
  end
  dflevel = figure('name','DFIELD6 Level sets',...
                 'vis','off',...
                 'numb','off','tag','DFIELD6');
  
  dfield6('figdefault',dflevel);
  hhsetup = get(0,'showhiddenhandles');
  set(dflevel,'menubar','none');
  
  lev.fr1 = uicontrol('style','frame');
  lev.fr2 = uicontrol('style','frame');
  
  inst1str = ['Enter the function in terms of the variables ',...
              tname, ' and ', xname,':'];

  lev.inst1 = uicontrol('style','text','horiz','left',...
                        'string',inst1str);
  
  lev.lfcn = uicontrol('style','edit','horiz','center',...
                       'string',lfcn,'call','',...
                       'background','w');
  
  lev.ch(3) = uicontrol('style','radio','horiz','center',...
                        'min',0,'max',3,...
                        'value',0,...
                        'vis','on',...
                        'string','Let DFIELD6 decide.');
  
  lev.ch(2) = uicontrol('style','radio','horiz','center',...
                        'min',0,'max',2,...
                        'value',0,...
                        'string','Select a point in the Display Window.');
      
  lev.inst2 = uicontrol('style','text','horiz','left',...
                        'string',['Choose one of the following ways to',...
                    ' choose level value(s):']);
      
  lev.ch(1)  = uicontrol('style','radio','horiz','center',...
                         'min',0,'max',1,...
                         'value',0,...
                         'string','Enter a vector of level values.');
      
  lev.rhs = uicontrol('style','edit','horiz','center',...
                      'string',' ','call','');

  lev.proc = uicontrol('style','push',...               
                          'string','Proceed',...
                          'call','dfield6(''levcomp'')');
  
  lev.close = uicontrol('style','push',...
                        'string','Close',...
                        'call','close');
      
  for i=1:3
    set(lev.ch(i),'user',lev.ch(:,[1:(i-1),(i+1):3]));
  end
  
  callrad = [
      'me = get(gcf,''currentobject'');',...
      'kk = get(me,''max'');',...
      'col = get(me,''backg'');',...
      'set(get(me,''user''),''value'',0),',...
      'set(me,''value'',kk);',...
      'ud = get(gcf,''user'');',...
      'if kk == 1,',...
      '   set(ud.rhs,''enable'',''on'',''backg'',''w'');',...
      'else,',...
      '   set(ud.rhs,''enable'',''off'',''backg'',col);',...
      'end,'];
  
  set(lev.ch,'call',callrad);
  
  left = 2; varl = 400; buttw = 60;
  nudge = 3;   
  tab = 15;
  lines1 = 5;
  lines2 = 2;
  xex = get(lev.inst1,'extent');
  ht = xex(4)+nudge;
  frw = varl + 2*tab+ 2*nudge;
  fr1bot = 2*left + ht;
  fr1ht = lines1*(nudge + ht) + nudge;
  fr2bot = fr1bot + fr1ht;
  fr2ht = lines2*(nudge + ht) + nudge;
  vbot = fr1bot + nudge;
  ch1bot = vbot + nudge + ht;
  ch2bot = ch1bot + nudge + ht;
  ch3bot = ch2bot + nudge + ht;
  inst2bot = ch3bot + nudge + ht;
  fbot =  fr2bot + nudge;
  inst1bot = fbot + nudge + ht;
  fleft = left + nudge + tab;
  instleft = left + nudge;
  chleft = instleft + tab;
  vleft = chleft + tab;
  vw = (frw - 2*vleft);
  fr1wind = [left,fr1bot,frw,fr1ht];
  fr2wind = [left,fr2bot,frw,fr2ht];
  inst1wind = [instleft,inst1bot,varl,ht];
  inst2wind = [instleft,inst2bot,varl,ht];
  fwind = [fleft,fbot,varl,ht];
  ch1wind = [chleft,ch1bot,varl,ht];
  ch2wind = [chleft,ch2bot,varl,ht];
  ch3wind = [chleft,ch3bot,varl,ht];
  vwind = [vleft,vbot,vw,ht];
  figw = 2*left + frw;
  fight = 3*left + ht + fr1ht + fr2ht;
  figwind = [40, 300, figw, fight];
  buttw = frw/2;
  sep = (figw - 2*buttw)/3;
  closel = sep;
  procl = 2*sep+buttw;
  clwind = [closel,left,buttw,ht];
  procwind = [procl,left,buttw,ht];
  set(dflevel,'pos',figwind);
  set(lev.fr1,'pos',fr1wind);
  set(lev.fr2,'pos',fr2wind);
  set(lev.inst1,'pos',inst1wind);
  set(lev.inst2,'pos',inst2wind);
  set(lev.ch(1),'pos',ch1wind);
  set(lev.ch(2),'pos',ch2wind);
  set(lev.ch(3),'pos',ch3wind);
  set(lev.rhs,'pos',vwind);
  set(lev.lfcn,'pos',fwind);
  set(lev.proc,'pos',procwind);
  set(lev.close,'pos',clwind);
  set(lev.ch(3),'value',3);
  set(lev.rhs,'enable','off');  
  
  child = get(dflevel,'children');
  set(dflevel,'vis','on','user',lev);
  set(child,'vis','on');
  
elseif strcmp(action,'levcomp')
  
  dflevel = gcf;
  ud = get(dflevel,'user');
  dfdisp = findobj('name','DFIELD6 Display');
  dud = get(dfdisp,'user');
  ch = ud.ch;
  val = zeros(1,3);
  for kk = 1:3
    val(kk) = get(ch(kk),'value');
  end
  KK = max(val);
  lfcn = get(ud.lfcn,'string');
  l=length(lfcn);
  for ( k = fliplr(findstr('.',lfcn)))
    if (find('*/^' == lfcn(k+1)))
      lfcn = [lfcn(1:k-1), lfcn(k+1:l)];
    end
    l=l-1;
  end
  l = length(lfcn);
  for (k=fliplr(find((lfcn=='^')|(lfcn=='*')|(lfcn=='/'))))
    lfcn = [lfcn(1:k-1) '.' lfcn(k:l)];
    l = l+1;
  end
  WINvect = dud.syst.wind;
  XxXxXx = WINvect(1) + rand*(WINvect(2)-WINvect(1));
  YyYyYy = WINvect(3) + rand*(WINvect(4)-WINvect(3));
  tname = dud.syst.tname;
  xname = dud.syst.xname;
  tname(find(abs(tname)==92))=[];  % Remove \s if any.
  xname(find(abs(xname)==92))=[];
  err = 0;res = 1;
  eval([tname,'=XxXxXx;'],'err = 1;');
  eval([xname,'=YyYyYy;'],'err = 1;');
  eval(['res = ',lfcn, ';'],'err = 1;');
  if err | isempty(res)
    errmsg = 'The function does not evaluate correctly.';
    fprintf('\a')
    errordlg(errmsg,'DFIELD6 error','on');
    return;
  end
  
  Xmin = WINvect(1);
  Xmax = WINvect(2);
  Ymin = WINvect(3);
  Ymax = WINvect(4);
  N = 50; k = 4;
  deltax=(Xmax - Xmin)/(N-1);
  deltay=(Ymax - Ymin)/(N-1);
  XXXg=Xmin + deltax*[-k:N+k];
  YYYg=Ymin + deltay*[-k:N+k];
  
  [Xx,Yy]=meshgrid(XXXg,YYYg);
  Xxx=Xx(:);Yyy=Yy(:);
  Ww = zeros(size(Xxx));
  eval([tname,'=Xxx;'],'err = 1;');
  eval([xname,'=Yyy;'],'err = 1;');
  eval(['Ww = ',lfcn, ';']);

  KKK = 3; %# of significant figures.
  
  switch KK
   case 1  % vector input
    rhs = get(ud.rhs,'string');
    rhs = str2num(rhs);
    
   case 2  % mouse input
    figure(dfdisp);
    [XX,YY] = ginput(1);
    figure(dflevel);
    eval([tname,'=XX;'],'err = 1;');
    eval([xname,'=YY;'],'err = 1;');
    eval(['rhs = ',lfcn, ';'],'err = 1'); 
    LL = ceil(log10(abs(rhs)));
    rhs = round(10^(KKK-LL)*rhs);
    rhs = 10^(LL-KKK)*rhs;
    
   case 3  % DFIELD6 input
    MM = max(Ww);mm = min(Ww);
    LL = ceil(log10(MM-mm));
    NN = 7;  % Number of curves
    rhs = mm+(1:NN).^2*(MM-mm)/NN^2;
    rhs = round(10^(KKK-LL)*rhs);
    rhs = 10^(LL-KKK)*rhs;
    
  end
  
  Ww = reshape(Ww,N+2*k+1,N+2*k+1);
  lrhs = length(rhs);
  if lrhs == 0
    return
  elseif lrhs == 1
    rhs = [rhs,rhs];
  end
  
  figure(dfdisp);
  [Cm,hcont] = contour(Xx,Yy,Ww,rhs,'--');
  hlabel = clabel(Cm,hcont);
%  set(hlabel,'fontsize',dud.fontsize,...
%            'color',dud.color.level,...
%            'rotation',0);
  set(hlabel,'fontsize',dud.fontsize,...
             'color',[1,0,0],...
             'rotation',0);
  set(hcont,'visible','on',...
            'color',dud.color.level,...
            'linestyle',':');
  dud.contours = [dud.contours ;hcont;hlabel];
  set(dfdisp,'user',dud);
  
elseif strcmp(action,'showbar')

  sbfig = gcbf;
  domymenu('menubar','toggletoolbar',sbfig);
  hhset = get(0,'showhiddenhandles');
  set(0,'showhiddenhandles','on');
  state = get(sbfig,'toolbar');
  if strcmp(state,'figure')
    fixtb = ['set(gcbo,''state'',''off'');'];
    set(findobj(sbfig,'tooltipstr','Zoom Out'),...
        'clickedcallback',['dfield6(''zoomback'');' fixtb]);
    set(findobj(sbfig,'tooltipstr','Zoom In'),...
        'clickedcallback',['dfield6(''zoomin'');' fixtb]);
    set(findobj(sbfig,'tooltipstr','Print'),...
        'clickedcallback','dfield6(''print'');');
    delete(findobj(sbfig,'tooltipstr','Rotate 3D'));
    sbmh = findobj(sbfig,'label','Show &Toolbar');
    set(sbmh,'label','Hide &Toolbar','checked','off');
  else
    sbmh = findobj(sbfig,'label','Hide &Toolbar');
    set(sbmh,'label','Show &Toolbar','checked','off');
  end
  set(0,'showhiddenhandles',hhset)
  
elseif strcmp(action,'figdefault')

  fig = input1;
  set(fig,'CloseRequestFcn','dfield6(''closefcn'')');
  dfset = findobj('name','DFIELD6 Setup');
  sud = get(dfset,'user');
  ud = get(fig,'user');
  ud.ssize = sud.ssize;
  fs = sud.fontsize;
  ud.fontsize = fs;
  style = sud.style;
  
  switch style
  case 'black'
    % if isunix | isvms, gamma = 0.5; else gamma = 0.0; end
    whitebg(fig,[0,0,0])
    if isunix | isvms
      fc = [.35 .35 .35];
    else
      fc = [.2 .2 .2];
    end
    set(fig,'color',fc);    
    set(fig,'defaultaxescolor',[0 0 0])
    % whitebg(fig,brighten([.2 .2 .2],gamma))
    set(fig,'defaultaxescolor',[0 0 0])
    set(fig,'defaultaxescolororder', ...
        1-[0 0 1;0 1 0;1 0 0;0 1 1;1 0 1;1 1 0;.25 .25 .25]) % ymcrgbw
    % set(fig,'colormap',brighten(jet(64),gamma))
    set(fig,'colormap',jet(64))
    set(fig,'defaultsurfaceedgecolor',[0 0 0]);
  case 'white'
    whitebg(fig,[1 1 1])
    set(fig,'color',[.8 .8 .8])
    set(fig,'defaultaxescolor',[1 1 1])
    set(fig,'defaultaxescolororder', ...
        [0 0 1;0 .5 0;1 0 0;0 .75 .75;.75 0 .75;.75 .75 0;.25 .25 .25]) % bgrymck
    set(fig,'colormap',jet(64))
    set(fig,'defaultsurfaceedgecolor',[0 0 0])

  case 'test'
    whitebg(fig,[1 1 1])
    set(fig,'color',[.8 .8 .8])
    set(fig,'defaultaxescolor',[1 1 1])
    set(fig,'defaultaxescolororder', ...
        [0 0 1;0 .5 0;1 0 0;0 .75 .75;.75 0 .75;.75 .75 0;.25 .25 .25]) % bgrymck
    set(fig,'colormap',jet(64))
    set(fig,'defaultsurfaceedgecolor',[0 0 0])
    % set(fig,'defaultuicontrolbackgroundcolor',[1 1 1]);

  case 'display'
    whitebg(fig,[1 1 1])
    set(fig,'defaultaxescolor',[1 1 1])
    set(fig,'defaultaxescolororder', ...
        [0 0 1;0 .5 0;1 0 0;0 .75 .75;.75 0 .75;.75 .75 0;.25 .25 .25]) % bgrymck
    set(fig,'colormap',jet(64))
    set(fig,'defaultsurfaceedgecolor',[0 0 0])
    set(fig,'color',[1 1 1]*220/255);
    set(fig,'defaultuicontrolbackgroundcolor',[1 1 1]*220/255);

  case 'bw'
    whitebg(fig,[0 0 0])
    set(fig,'color',[0 0 0])
    set(fig,'defaultaxescolor',[0 0 0])
    set(fig,'defaultaxescolororder', ...
        [1 1 1])
    set(fig,'colormap',[1 1 1;0 0 0])
    set(fig,'defaultsurfaceedgecolor',[0 0 0])
    
  end 
  set(fig,'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);
  
  
elseif strcmp(action,'paraeval')
  
  %     Replace a parameter with its value in string form.
  
  para = deblank(input1);
  value = input2;
  value = ['(',value,')'];
  str = input3;
  ll = length(str);
  lp = length(para);
  lv = length(value);
  
  if strcmp(para,str)
    str = value
  elseif (ll >= lp+1)
    k = findstr(para,str);
    
    lk = length(k);
    lopstr = '(+-*/^';
    ropstr = ')+-*/^';
    s = [];
    pos = 1;
    for jj = 1:lk
      if (((k(jj) == 1)|(find(lopstr == str(k(jj)-1))))...
            &((k(jj)+lp-1 == ll)|(find(ropstr == str(k(jj) + lp)))))
        s = [s,str(pos:(k(jj)-1)),value];
        pos = k(jj)+lp;
      end
    end
    str = [s,str(pos:ll)];
  end
  output = str;

elseif strcmp(action,'quit')
  
  dfset = findobj('name','DFIELD6 Setup');
  sud = get(dfset,'user');
  if sud.remtd 
    rmpath(tempdir);
  end
  oldfiles = dir([tempdir,'dftp*.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 DFIELD6')
      delete(fn)
    end
  end
  h = findobj('tag','DFIELD6');
  delete(h);
  
elseif strcmp(action,'restart')
  
  dfdisp = findobj('name','DFIELD6 Display');
  dfset = findobj('name','DFIELD6 Setup');
  if (~isempty(dfdisp))
    dud = get(dfdisp,'user');
    if isfield(dud,'function')
      fcn = dud.function;
      if (exist(fcn)==2) delete([fcn,'.m']);end
    end
  end
  h = findobj('tag','DFIELD6');
  delete(setdiff(h,dfset));
  sud = get(dfset,'user');
  sud.flag = 0;
  set(dfset,'user',sud);
  figure(dfset)
  
elseif strcmp(action,'closefcn')
  
  fig = gcf;  
  name = get(fig,'name');
  if strcmp(name,'DFIELD6 Setup') | strcmp(name,'DFIELD6 Display')
    quest = ['Closing this window will cause all DFIELD6 ',...
          'windows to close, and DFIELD6 will stop.  ',...
          'Do you want to quit DFIELD6?'];
    butt = questdlg(quest,'Quit DFIELD6?','Quit','Cancel','Quit');
    if strcmp(butt,'Quit')
      dfield6('quit');
    end
  elseif strcmp(name,'DFIELD6 Linearization')
    dud = get(fig,'user');
    fcn = dud.function;
    if (exist(fcn)==2) delete([fcn,'.m']);end    
    delete(findobj('label',name));
    delete(fig);
  else
    delete(findobj('label',name));
    delete(fig);
  end
  
elseif strcmp(action,'text')
   
   dfdisp = gcf;
   oldcall = get(dfdisp,'WindowButtonDownFcn');
   set(dfdisp,'WindowButtonDownFcn','');
   prompt = ['Enter the text here. Then choose ',...
         'the location in the Display Window.'];
   txtstr = inputdlg(prompt,'Text entry');
   if ~isempty(txtstr)
      txtstr = txtstr{1};
      figure(dfdisp);
      gtext(txtstr);
   end
   set(dfdisp,'WindowButtonDownFcn',oldcall);
   
elseif strcmp(action,'confused')
  
  tstring = 'DFIELD6 is totally confused';
  qstring = {['You will have to restart DFIELD6 from '...
        'the beginning in order to ',...
        'do anything new.  However, it might be possible '...
        'to save the current system ',...
        'or the gallery to make your restart easier, '...
        'or it may be possible to ',...
        'print out a figure, if the appropriate '...
        'figures are visible.  In such a case ',...
        'it would be best to do nothing now.'];
        'What do you want to do?'};
  bstr1 = 'Quit and restart DFIELD6.';
  bstr2 = 'Just quit DFIELD6.';
  bstr3 = 'Do nothing.';
  answer = questdlg(qstring,tstring,bstr1,bstr2,bstr3,bstr1);
  if strcmp(answer,bstr1)
    delete(findobj('tag','DFIELD6'));
    dfield6;return
  elseif strcmp(answer,bstr2)
    delete(findobj('tag','DFIELD6'));
    return
  else
    return
  end
  
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [tout,yout] = dfdp45(dfcn,tspan,y0,disph)

% DFDP45 is an implementation of the explicit Runge-Kutta (4,5) which
%        is described in Chapters 5 & 6 of John Dormand's book, 
%        Numerical Methods for Differential Equations.  
% 
%        This is the same algorithm used in ODE45, part of the new MATLAB 
%        ODE Suite.  Details are to be found in The MATLAB ODE Suite, 
%        L. F. Shampine and M. W. Reichelt, SIAM Journal on Scientific 
%        Computing, 18-1, 1997.


% Input the user data. 

dud = get(disph,'user');
dispha = dud.axes;
ud = get(dispha,'user');
refine = dud.settings.refine;
tol = dud.settings.tol;
hmax = dud.settings.hmax;
DY = ud.DY;
col = dud.color.temp;
speed = dud.settings.speed;
slow = (speed < 100);


% Initialize the stopping criteria.


%  The compute window.

CC = ud.cwind(3);
DD = ud.cwind(4);

% The stop button.

stop = 0;
ud.stop = 0;

% Set the the line handle.

ph = plot([tspan(1),tspan(1)],[y0,y0],'color',col,...
    'erase','none',...
    'parent',dispha);
ud.line = ph;

% Set up the phase line

vi = dud.pline;
figure(disph);
v = axis;
aa = v(1)+0.01*(v(2)-v(1));
plh = plot(aa,y0,'.','markersize',20,'color',col,...
    'erase','xor','parent',dispha,'visible',vi);    
ud.pline = plh;
set(dispha,'UserData',ud);

% Initialize the loop.

t0 = tspan(1);
tfinal = tspan(2);
tdir = sign(tfinal - t0);
t = t0;
y = y0;

% By default, hmax is 1/10 of the interval.
%hmax = min(abs(0.1*(tfinal-t)),1);

rDY = DY(:,ones(1,refine));
steps = 0;
block = 120;
tout = zeros(block,1);
yout = zeros(block,1);

N = 1;
tout(N) = t;
yout(N) = y;

% Initialize method parameters.
pow = 1/5;
C = [1/5; 3/10; 4/5; 8/9; 1; 1]; 
A = [
    1/5         3/40    44/45   19372/6561      9017/3168  
    0           9/40    -56/15  -25360/2187     -355/33    
    0           0       32/9    64448/6561      46732/5247      
    0           0       0       -212/729        49/176          
    0           0       0       0               -5103/18656     
    0           0       0       0               0               
    0           0       0       0               0            
    ];
bhat = [35/384 0 500/1113 125/192 -2187/6784 11/84 0]';
% E = bhat - b.
E = [71/57600; 0; -71/16695; 71/1920; -17253/339200; 22/525; -1/40];
if refine > 1
  sigma = (1:refine-1) / refine;
  S = cumprod(sigma(ones(4,1),:));
  bstar = [
    1       -183/64     37/12       -145/128
    0       0       0       0
    0       1500/371    -1000/159   1000/371
    0       -125/32     125/12      -375/64 
    0       9477/3392   -729/106    25515/6784
    0       -11/7       11/3        -55/28
    0       3/2     -4      5/2
    ];
  bstar = bstar*S;
    
end

f = zeros(1,7);
f0 = feval(dfcn,t,y);

mm = max(abs(f0./DY));
absh = hmax;
if mm
  absh = min(absh,1/(100*mm));
end

f(:,1) = f0;


% THE MAIN LOOP
tic

while ~stop
  
  % hmin is a small number such that t+h is distinquishably
  % different from t if abs(h) > hmin.
  hmin = 16*eps*abs(t);
  absh = min(hmax, max(hmin, absh));    
  h = tdir * absh;
    
  % LOOP FOR ADVANCING ONE STEP.
  while stop~=5
    hC= h * C;
    hA = h * A;

    f(:,2) = feval(dfcn, t + hC(1), y + f*hA(:,1));
    f(:,3) = feval(dfcn, t + hC(2), y + f*hA(:,2));
    f(:,4) = feval(dfcn, t + hC(3), y + f*hA(:,3));
    f(:,5) = feval(dfcn, t + hC(4), y + f*hA(:,4));
    f(:,6) = feval(dfcn, t + hC(5), y + f*hA(:,5));
    tn = t + h;
    yn = y + f*h*bhat;
    f(:,7) = feval(dfcn, tn, yn);
    
    % Estimate the error.
    err = abs(h * f * E);
    alpha = (2*max(err/((abs(y)+abs(yn)+DY(2)*1e-3)*tol)))^pow;
    if alpha < 1         % Step is OK
      break
    else
      if absh <= hmin           % This keeps us out of an infinite loop.
        stop = 5;
        break;
      end
      
      absh = max(hmin,0.8*absh / alpha);
      h = tdir * absh;
    end  % if alpha < 1
  end  % while stop~=5
  steps = steps + 1;
  

  oldN = N;
  N = N + refine;
  if N > length(tout)
    tout = [tout; zeros(block,1)];  
    yout = [yout; zeros(block,1)];
  end
  if refine > 1             % computed points, with refinement
    j = oldN+1:N-1;
    tout(j) = t + h*sigma';
    yout(j) = (y(:,ones(length(j),1)) + f*h*bstar).';
    tout(N) = tn;
    yout(N) = yn;
  else               % computed points, no refinement
    tout(N) = tn;
    yout(N) = yn;
  end
  
  % Update stop.   Maybe the stop button has been pressed.

  ud = get(dispha,'user');
  stop = max(ud.stop,stop);
  
  % Are we out of the compute window?
  
  if yn<CC | yn>DD
    stop = 1;
  end  
  if (abs(tn-tfinal) < hmin)
    stop = 2;
  end  % if gstop
  
  i = oldN:N;    
  set(ph,'Xdata',tout(i),'Ydata',yout(i));
  set(plh,'Xdata',aa,'Ydata',yn);
  drawnow    
  
  % Compute a new step size.
  absh = max(hmin,0.8*absh / max( alpha,0.1));
  absh = min(absh,tdir*(tfinal - tn));
  h = absh*tdir;
  % Advance the integration one step.
  t = tn;
  y = yn;
  f(1) = f(7);                      % Already evaluated dfcn(tnew,ynew)
  
  if slow
    ttt= N/(speed*refine);
    tt = toc;
    while tt < ttt
      tt = toc;
    end
  end
end  % while ~stop
tout = tout(1:N);
yout = yout(1:N,:);

if dud.notice
  nstr = get(dud.notice,'string');
  
  switch stop
    case 1
      nstr{5} = [nstr{5}, ' left the computation window.']; 
    case 4
      nstr{5} = [nstr{5}, ' was stopped by the user.'];
    case 5 
      ystr = ['(',num2str(t,2), ', ', num2str(y,2),').'];
      nstr(1:3) = nstr(2:4);
      nstr{4} = [nstr{5},' experienced a failure at ',ystr];
      nstr{5} = 'Problem is singular or tolerances are too large.';
  end
  set(dud.notice,'string',nstr);
end
drawnow

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [tout,yout] = dfeul(dfcn,tspan,y0,disph)

% DFEUL  is an implementation of Euler's method


% Input the user data. 

dud = get(disph,'user');
dispha = dud.axes;
ud = get(dispha,'user');
refine = dud.settings.refine;
ssize = dud.settings.stepsize;
DY = ud.DY;
col = dud.color.temp;
speed = dud.settings.speed;
slow = (speed < 100);

% Initialize the stopping criteria.

%  The compute window.

CC = ud.cwind(3);
DD = ud.cwind(4);

% The stop button.

stop = 0;
ud.stop = 0;

% Set the the line handle.

ph = plot([tspan(1),tspan(1)],[y0,y0],'color',col,...
    'erase','none',...
    'parent',dispha);
ud.line = ph;

% Set up the phase line

vi = dud.pline;
figure(disph);
v = axis;
aa = v(1)+0.01*(v(2)-v(1));
plh = plot(aa,y0,'.','markersize',20,'color',col,...
    'erase','xor','parent',dispha,'visible',vi);    
ud.pline = plh;
set(dispha,'UserData',ud);

% Initialize the loop.

t0 = tspan(1);
tfinal = tspan(2);
tdir = sign(tfinal - t0);
t = t0;
y = y0;

h = ssize*tdir;
block = 120;
tout = zeros(block,1);
yout = tout;
N = 1;
tout(N) = t;
yout(N) = y;

% The main loop
tic
while ~stop
  if abs(t - tfinal) < ssize 
    h = tfinal - t; 
  end

  % Compute the slope
  s1 = feval(dfcn,t,y); 
  t = t + h;
  y = y + h*s1;         
  if N >= length(tout)
    tout = [tout;zeros(block,1)];
    yout = [yout;zeros(block,1)];
  end
  N = N + 1;  
  tout(N) = t;
  yout(N) = y;
  
  % Update stop.   Maybe the stop button has been pressed.

  ud = get(dispha,'user');
  stop = max(ud.stop,stop);
  
  % Are we out of the compute window?
  
  if y<CC | y>DD
    stop = 1;
  end  
  if (abs(t-tfinal) < 0.0001*ssize)
    stop = 2;
  end  
  
  i = (N-1):N;    
  set(ph,'Xdata',tout(i),'Ydata',yout(i));
  set(plh,'Xdata',aa,'Ydata',y);
  drawnow    
  if slow
    ttt= N/(speed);
    tt = toc;
    while tt < ttt
      tt = toc;
    end
  end
    
end  % while ~stop

tout = tout(1:N);
yout = yout(1:N,:);
if dud.notice
  nstr = get(dud.notice,'string');
  
  switch stop
    case 1
      nstr{5} = [nstr{5}, ' left the computation window.']; 
    case 4
      nstr{5} = [nstr{5}, ' was stopped by the user.'];
    case 5 
      ystr = ['(',num2str(t,2), ', ', num2str(y,2),').'];
      nstr(1:3) = nstr(2:4);
      nstr{4} = [nstr{5},' experienced a failure at ',ystr];
      nstr{5} = 'Problem is singular or tolerances are too large.';
  end
  set(dud.notice,'string',nstr);
end
drawnow


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [tout,yout] = dfrk2(dfcn,tspan,y0,disph)

% DFRK2  is an implementation of the second order Runge-Kutta method. 


% Input the user data. 

dud = get(disph,'user');
dispha = dud.axes;
ud = get(dispha,'user');
refine = dud.settings.refine;
ssize = dud.settings.stepsize;
DY = ud.DY;
col = dud.color.temp;
speed = dud.settings.speed;
slow = (speed < 100);

% Initialize the stopping criteria.

%  The compute window.

CC = ud.cwind(3);
DD = ud.cwind(4);

% The stop button.

stop = 0;
ud.stop = 0;

% Set the the line handle.

ph = plot([tspan(1),tspan(1)],[y0,y0],'color',col,...
    'erase','none',...
    'parent',dispha);
ud.line = ph;

% Set up the phase line

vi = dud.pline;
figure(disph);
v = axis;
aa = v(1)+0.01*(v(2)-v(1));
plh = plot(aa,y0,'.','markersize',20,'color',col,...
    'erase','xor','parent',dispha,'visible',vi);    
ud.pline = plh;
set(dispha,'UserData',ud);

% Initialize the loop.

t0 = tspan(1);
tfinal = tspan(2);
tdir = sign(tfinal - t0);
t = t0;
y = y0;

h = ssize*tdir;
block = 120;
tout = zeros(block,1);
yout = tout;
N = 1;
tout(N) = t;
yout(N) = y;

% The main loop
tic
while ~stop
  if abs(t - tfinal) < ssize 
    h = tfinal - t; 
  end

  % Compute the slope
  s1 = feval(dfcn,t,y); 
  s2 = feval(dfcn, t + h, y + h*s1); 
  t = t + h;
  y = y + h*(s1 + s2)/2;

  if N >= length(tout)
    tout = [tout;zeros(block,1)];
    yout = [yout;zeros(block,1)];
  end
  N = N + 1;  
  tout(N) = t;
  yout(N) = y;
  
  % Update stop.   Maybe the stop button has been pressed.

  ud = get(dispha,'user');
  stop = max(ud.stop,stop);
  
  % Are we out of the compute window?
  
  if y<CC | y>DD
    stop = 1;
  end  
  if (abs(t-tfinal) < 0.0001*ssize)
    stop = 2;
  end  
  
  i = (N-1):N;    
  set(ph,'Xdata',tout(i),'Ydata',yout(i));
  set(plh,'Xdata',aa,'Ydata',y);
  drawnow
  if slow
    ttt= N/(speed);
    tt = toc;
    while tt < ttt
      tt = toc;
    end
  end
    
end  % while ~stop

tout = tout(1:N);
yout = yout(1:N,:);
if dud.notice
  nstr = get(dud.notice,'string');
  
  switch stop
    case 1
      nstr{5} = [nstr{5}, ' left the computation window.']; 
    case 4
      nstr{5} = [nstr{5}, ' was stopped by the user.'];
    case 5 
      ystr = ['(',num2str(t,2), ', ', num2str(y,2),').'];
      nstr(1:3) = nstr(2:4);
      nstr{4} = [nstr{5},' experienced a failure at ',ystr];
      nstr{5} = 'Problem is singular or tolerances are too large.';
  end
  set(dud.notice,'string',nstr);
end
drawnow



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [tout,yout] = dfrk4(dfcn,tspan,y0,disph)

% DFRK4  is an implementation of the fourth order Runge-Kutta method. 


% Input the user data. 

dud = get(disph,'user');
dispha = dud.axes;
ud = get(dispha,'user');
refine = dud.settings.refine;
ssize = dud.settings.stepsize;
DY = ud.DY;
col = dud.color.temp;
speed = dud.settings.speed;
slow = (speed < 100);

% Initialize the stopping criteria.

%  The compute window.

CC = ud.cwind(3);
DD = ud.cwind(4);

% The stop button.

stop = 0;
ud.stop = 0;

% Set the the line handle.

ph = plot([tspan(1),tspan(1)],[y0,y0],'color',col,...
    'erase','none',...
    'parent',dispha);
ud.line = ph;

% Set up the phase line

vi = dud.pline;
figure(disph);
v = axis;
aa = v(1)+0.01*(v(2)-v(1));
plh = plot(aa,y0,'.','markersize',20,'color',col,...
    'erase','xor','parent',dispha,'visible',vi);    
ud.pline = plh;
set(dispha,'UserData',ud);

% Initialize the loop.

t0 = tspan(1);
tfinal = tspan(2);
tdir = sign(tfinal - t0);
t = t0;
y = y0;

h = ssize*tdir;
block = 120;
tout = zeros(block,1);
yout = tout;
N = 1;
tout(N) = t;
yout(N) = y;

% The main loop
tic
while ~stop
  if abs(t - tfinal) < ssize 
    h = tfinal - t; 
  end

  % Compute the slope
  s1 = feval(dfcn,t,y); 
  s2 = feval(dfcn, t + h/2, y + h*s1/2); s2=s2(:);
  s3 = feval(dfcn, t + h/2, y + h*s2/2); s3=s3(:);
  s4 = feval(dfcn, t + h, y + h*s3); s4=s4(:);
  
  t = t + h;
  y = y + h*(s1 + 2*s2 + 2*s3 +s4)/6;

  if N >= length(tout)
    tout = [tout;zeros(block,1)];
    yout = [yout;zeros(block,1)];
  end
  N = N + 1;  
  tout(N) = t;
  yout(N) = y;
  
  % Update stop.   Maybe the stop button has been pressed.

  ud = get(dispha,'user');
  stop = max(ud.stop,stop);
  
  % Are we out of the compute window?
  
  if y<CC | y>DD
    stop = 1;
  end  
  if (abs(t-tfinal) < 0.0001*ssize)
    stop = 2;
  end  
  
  i = (N-1):N;    
  set(ph,'Xdata',tout(i),'Ydata',yout(i));
  set(plh,'Xdata',aa,'Ydata',y);
  drawnow    
  if slow
    ttt= N/(speed);
    tt = toc;
    while tt < ttt
      tt = toc;
    end
  end
    
end  % while ~stop

tout = tout(1:N);
yout = yout(1:N,:);
if dud.notice
  nstr = get(dud.notice,'string');
  
  switch stop
    case 1
      nstr{5} = [nstr{5}, ' left the computation window.']; 
    case 4
      nstr{5} = [nstr{5}, ' was stopped by the user.'];
    case 5 
      ystr = ['(',num2str(t,2), ', ', num2str(y,2),').'];
      nstr(1:3) = nstr(2:4);
      nstr{4} = [nstr{5},' experienced a failure at ',ystr];
      nstr{5} = 'Problem is singular or tolerances are too large.';
  end
  set(dud.notice,'string',nstr);
end
drawnow

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%