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

% DFIELD5	is an interactive tool for studying single first order
%	differential equations.  When DFIELD5 is executed, a DFIELD5 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 DFIELD5 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) John C. Polking, Rice University

%  Version 5.01
%  Last modified:  September 4, 1997

style = 'white';
startstr = 'dfield5';
ud.fontsize = get(0,'defaultaxesfontsize');

if nargin < 1
  action ='initialize';
end


%if isa(action,'numeric')
%  ud.fontsize = max(5,min(action,25));
  % set(0,'defaultaxesfontsize',action);
%  action   ='initialize';
%end  

% disp(action)
if strcmp(action,'black')
  action   ='initialize';
  style = 'black';
  startstr = 'dfield5(''black'')';
elseif strcmp(action,'display')
  action   ='initialize';
  style = 'display';
  startstr = 'dfield5(''display'')';
end

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

if strcmp(action,'initialize')

  % First we make sure that there is no other copy of DFIELD5
  % running, since this causes problems.	
  dfh = findobj('tag','DFIELD5');
  if ~isempty(dfh);
    qstring = {'There are some open DFIELD5 figures (perhaps invisible).  ';...
    'What do you want to do?'};
    tstring = 'Only one copy of DFIELD5 can be open at one time.';
    b1str = 'Restart DFIELD5.';
    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);
  
  % Next we look for old files created by DFIELD5.
  
  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 DFIELD5')
      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

  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 = 20;  
  ud.style = style;
  ud.settings.magn = 1.25;
  ud.settings.refine = 4;
  ud.settings.tol = 5e-4;
  ud.settings.solver = 'Dormand-Prince';
  ud.settings.stepsize = 0.1;
  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 = [0 1 1]; % cyan for arrows
    color.arrows = .5*[1 1 1];  % gray for arrows
    color.pline = [1 1 1];    
   
  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];    
    
  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];    
    
  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];    
  end
  ud.color = color;
  
  dfset = figure('name','DFIELD5 Setup','numb','off',...
      'tag','DFIELD5','visible','off',...
      'user',ud);
  
  dfield5('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.
  if strcmp(style,'display')
    dfigurebot = 230;
  else
    dfigurebot = 30;	% Bottom of the figure.
  end
  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]);
  
  ud.h.xname=uicontrol('pos',[eqleft, derbot, varw, texth],'style','edit',...
      'horizon','right','string',ud.o.xname,'call',xname,'visible','off');
  
  eq(2) = uicontrol('style','text',...
      'pos',[eqleft+varw derbot-fudge equalw texth],...
      'horizon','center','string',''' = ','visible','off');
  
  ud.h.der=uicontrol('pos',[eqleft+varw+equalw derbot eqlength texth],...
      'string',ud.o.der,...
      'horizon','left','style','edit',...
      'call',der,'visible','off');
  
  eq(3) = uicontrol('style','text',...
      'horizon','right','string','The independent variable is ',...
      'visible','off');
 
  ext = get(eq(3),'extent');
  tnw = ext(3)+2;
  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');
  
  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 = 40*rr;
  pvalw = 50*rr;
  peqw = 10*rr;
  pbot = pfrbot + 5;

  pleft1 = eqleft + 50*rr + 5;
  peqleft1 = pleft1 + pnamew;
  pvleft1 = peqleft1 + peqw;
  pleft2 = dfigwidth - 10 - pnamew - pvalw - peqw;
  peqleft2 = pleft2 + pnamew;
  pvleft2 = peqleft2 + peqw;	
  
  paratit=uicontrol('style','text',...
      'pos',[eqleft pbot-fudge 50*rr texth],...
      'horizon','center',...
      'string','Parameters:','visible','off');
  
  ud.h.pname1 = uicontrol('style','edit',...
      'pos',[pleft1 pbot pnamew texth],...
      'horizon','right','string',ud.o.pname1,...
      'call',p1,'visible','off');
  
  equal(1) = uicontrol('style','text',...
      'pos',[peqleft1 pbot-fudge peqw texth],...
      'horizon','center','string','=','visible','off');
  
  ud.h.pval1 = uicontrol('style','edit',...
      'pos',[pvleft1 pbot pvalw texth],...
      'string',ud.o.pval1,'call',p1,'visible','off');
  
  ud.h.pname2 = uicontrol('style','edit',...
      'pos',[pleft2 pbot pnamew texth],...
      'horizon','right','string',ud.o.pname2,...
      'call',p2,'visible','off');
  
  equal(2) = uicontrol('style','text',...
      'pos',[peqleft2 pbot-fudge peqw texth],...
      'horizon','center','string','=','visible','off');
  
  ud.h.pval2 = uicontrol('style','edit',...
      'pos',[pvleft2 pbot pvalw texth],...
      'string',ud.o.pval2,'call',p2,'visible','off');
  
  
  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');
  
  % 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');
  
  ud.h.wind(1) = uicontrol('style','edit',...
      'pos',[eqleft+twindw winbot2 windw texth],...
      'string',num2str(ud.o.wind(1)),...
      'call',w1,'visible','off');
  
  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');
  
  ud.h.wind(2) = uicontrol('style','edit',...
      'pos',[eqleft+twindw winbot3 windw texth],...
      'string',num2str(ud.o.wind(2)),...
      'call',w2,'visible','off');
  
  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');
  
  ud.h.wind(3) = uicontrol('style','edit',...
      'pos',[twindl+twindw winbot2 windw texth],...
      'string',num2str(ud.o.wind(3)),...
      'call',w3,'visible','off');
  
  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');
  
  ud.h.wind(4) = uicontrol('style','edit',...
      'pos',[twindl+twindw winbot3 windw texth],...
      'string',num2str(ud.o.wind(4)),...
      'call',w4,'visible','off');
  
  
  
  butt(1) = uicontrol('style','push',...
      'pos',qwind,...
      'string','Quit','call','dfield5(''quit'')','visible','off');
  
  butt(2) = uicontrol('style','push',...
      'pos',rwind,...
      'string','Revert','call','dfield5(''revert'')','visible','off');
  
  butt(3) = uicontrol('style','push',...
      'pos',pwind,...
      'string','Proceed','call','dfield5(''proceed'')','visible','off');
  
  
  mefile = findobj(dfset,'label','File');

  meedit = findobj(gcf,'label','Edit'); % uimenu('label','Edit');
  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');
  merevert = uimenu(mefile,'label','Revert','call',...
      'dfield5(''revert'')','separator','on');
  meproceed = uimenu(mefile,'label','Proceed','call',...
      'dfield5(''proceed'')','separator','on',...
      'accelerator','G');
  merestart = uimenu(mefile,'label','Restart DFIELD5','call',...
      'dfield5(''restart'')','separator','on');
  mequit = uimenu(mefile,'label','Quit DFIELD5','call',...
      'dfield5(''quit'')','separator','on');
  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,'DFIELD5 error','on');
    return
  end
  xstr = sud.c.xname;
  if isempty(xstr)
    errmsg = 'The first dependent variable needs a name.';
    fprintf('\a')
    errordlg(errmsg,'DFIELD5 error','on');
    return
  end
  tstr = sud.c.tname;
  if isempty(tstr)
    errmsg = 'The second dependent variable needs a name.';
    fprintf('\a')
    errordlg(errmsg,'DFIELD5 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,'DFIELD5 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,'DFIELD5 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','DFIELD5 Display');
    dud = get(dfdisp,'user');
    aud = get(dud.axes,'user');
    tstr = get(get(dud.axes,'title'),'string');
    wind = sud.c.wind(:);
    if (~all(wind == dud.syst.wind(:)))
      dwind = [wind(1); wind(3); -wind(2); -wind(4)];
      DY = [wind(2)-wind(1); wind(4)-wind(3)];
      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;
    set(dfdisp,'user',dud);
    dfield5('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))=[];
    
    % 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
        
    txderstr = derivstr;
  
    kxder = find(abs(txderstr)==42);
    % txderstr = derivstr;
    txderstr(kxder)=' '*ones(size(kxder));
    txderstr = strrep(txderstr,'-',' - ');
    txderstr = strrep(txderstr,'+',' + ');
    if (abs(txderstr(1)) == 32)
      txderstr = txderstr(2:length(txderstr));
    end
    tstr = [Xname,' '' = ', txderstr];
    
    % Replace the parameters with their values.
    pflag = zeros(1,2);
    perr = [];
    if ~isempty(pname1)
      if isempty(str2num(parav1))
	perr = sud.h.pval1;
      else
	xderivstr = dfield5('paraeval',pname1,parav1,derivstr);
	if ~strcmp(xderivstr,derivstr)
	  pflag(1) = 1;
	  derivstr = xderivstr;
	end
      end
    end    
    if ~isempty(pname2)
      if isempty(str2num(parav2))
	perr = [perr;sud.h.pval2]
      else
	xderivstr = dfield5('paraeval',pname2,parav2,derivstr);
	if ~strcmp(xderivstr,derivstr)
	  pflag(2) = 1;
	  derivstr = xderivstr;
	end
      end      
    end

    % 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))=[];
    
    SS = warning;
    warning off
    XxXxXx = WINvect(1) + rand*(WINvect(2)-WINvect(1));
    TtTtTt = WINvect(3) + rand*(WINvect(4)-WINvect(3));
    err = 0;
    eval([Xname,'=XxXxXx;'],'err = 1;');
    if err
      errmsg = ['"',xstr, '" is not a valid variable name in MATLAB.'];
      fprintf('\a')
      errordlg(errmsg,'DFIELD5 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,'DFIELD5 error','on');
      warning(SS)
      return
    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
    
    err = 0;res = 1;
    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,'DFIELD5 error','on');
      dfield5('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 DFIELD5\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([dfcn,'.m'],'w');
    fprintf(dff,fcnstr);
    fprintf(dff,commstr);
    fprintf(dff,varstr);
    fprintf(dff,lenstr);
    fprintf(dff,derstr1);
    fprintf(dff,derstr2);
    fclose(dff);


    % Find DFIELD5 Display if it exists.	
    % If DFIELD5 Display exists, update it.  If it does not build it.
    
    dfdisp = findobj('name','DFIELD5 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
      dud.syst = sud.c;
      dud.settings = sud.settings;
      dfdisp = figure('name','DFIELD5 Display',...
	  'numb','off','interrupt','on',...
	  'visible','off','tag','DFIELD5');
      dfield5('figdefault',dfdisp);
      fs = get(dfdisp,'defaultaxesfontsize');
      r = fs/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
      nframeh = 70;   % Default notice frame height
      dfaxb = 4+nframeh+35;

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

      uni = get(0,'units');
      set(0,'units','pixels');
      ss = get(0,'screensize');
      set(0,'units',uni);
      sw = ss(3);sh = ss(4);
      if strcmp(sud.style, 'display')
	bottom = 230;
      else 
	bottom = 10;
      end
      
      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+10*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 - 35;
      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','dfield5(''quit'')','visible','off');
      
      dbutt(3) = uicontrol('style','push',...
	  'pos',[buttl,dfaxb+buttsep,buttw,butth],...
	  'string','Print',...
	  'call','dfield5(''print'')','visible','off');	
      
      dud.butt = dbutt;
      
      fmenu = findobj(gcf,'label','File');  % uimenu('label','File');
      emenu = findobj(gcf,'label','Edit'); % uimenu('label','Edit');
      
      menu(1) = uimenu('label','Options','visible','off');
      
      meprint = uimenu(fmenu,'label','Print DFIELD5 Display',...
	  'Handlevis','on', ...
	  'Accelerator','P', ...
	  'Callback','print -noui');
      
      meprsetup =   uimenu(fmenu,'label','Page Position ...',...
	  'Handlevis','on', ...
	  'Callback','pagedlg(gcf)'); 
      
      meprdlg =   uimenu(fmenu,'label','Print ...',...
	  'Handlevis','on', ...
	  'Callback','printdlg(gcf)');  
      
      mequit = uimenu(fmenu,'label','Quit DFIELD5',...
	  'call','dfield5(''quit'')','separator','on');
      
      menukey = uimenu(menu(1),'label','Keyboard input.','call',...
	  'dfield5(''kbd'')');
      
      mesev   = uimenu(menu(1),'label','Plot several solutions.',...
	  'call','dfield5(''several'')');
      
      menu(2) = uimenu(emenu,'label','Zoom in.',...
	  'call','dfield5(''zoomin'')');
      
      zbmenu = uimenu(emenu,'label','Zoom back.','call',...
	  'dfield5(''zoomback'')','enable','off','tag','zbmenu');
      
      medall = uimenu(emenu,'label','Erase all solutions.',...
	  'call','dfield5(''dall'')','separator','on');
      
      medel = uimenu(emenu,'label','Delete a graphics object.',...
	  'call','dfield5(''delete'')','visible','on');
      
      menutext = uimenu(emenu,'label','Enter text on the Display Window.',...
	  'call','dfield5(''text'')');
      
      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'');',...	  
	    '  if ~isempty(ud.plhand);',...	  
	    '    delete(ud.plhand);',...
	    '    ud.plhand = [];',...	  
	    '  end,',...	  
	    '  ud.pland =[];',...
	    '  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','dfield5(''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;',...
	    'set(ud.solver,''checked'',''off'');',...
	    'set(me,''checked'',''on'');',...
	    'set(gcf,''user'',ud);',...
	    'dfset = findobj(''name'',''DFIELD5 Setup'');',...
	    'sud = get(dfset,''user'');',...
	    'sud.settings = ud.settings;',...
	    'set(dfset,''user'',sud);',...
	    'dfield5(''solvset'');'];
      dpset.refine = 4;
      dpset.tol = 5e-4;
      dpset.solver = 'Dormand-Prince';
      dpset.stepsize = 0;
      eulset.refine = 1;
      eulset.tol = 0;
      eulset.solver = 'Euler';
      eulset.stepsize = 0.1;
      rk2set.refine = 1;
      rk2set.tol = 0;
      rk2set.solver = 'Runge-Kutta 2';
      rk2set.stepsize = 0.1;
      rk4set.refine = 1;
      rk4set.tol = 0;
      rk4set.solver = 'Runge-Kutta 4';
      rk4set.stepsize = 0.1;
      
      dud.solver(1) = uimenu(mesolve,'label','Dormand-Prince',...
	  'checked','on',...
	  'call',solverstr,'user',dpset);
      
      dud.solver(2) = uimenu(mesolve,'label','Euler',...
	  'checked','off',...
	  'call',solverstr,'user',eulset);
      
      dud.solver(3) = uimenu(mesolve,'label','Runge-Kutta 2',...
	  'checked','off',...
	  'call',solverstr,'user',rk2set);
      
      dud.solver(4) = uimenu(mesolve,'label','Runge-Kutta 4',...
	  'checked','off',...
	  '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','dfield5(''settings'')');
      
      menu(6) = uimenu(menu(1),'label','Make the Display Window inactive.',...
	  'call','dfield5(''hotcold'')','separator','on');
      
      dud.menu = menu;
      
      if exist('RICE') == 2
	ricemenu = uimenu (menu(1),'label','Set default printer.',...
	    'call','dfield5 (''riceprint'')','separator','on');
      end
      
      
      set(gcf,'WindowButtonDownFcn','dfield5(''down'')');
      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 = [];    
    end   % if (~isempty(dfdisp))  & else
    
    dfdispa = dud.axes;
    axes(dfdispa);
    xlabel(Tname);
    ylabel(Xname);

    set(findobj('tag','zbmenu'),'enable','off');
    k = find(pflag);
    if ~isempty(k)
      lk = length(k);
      pstring = cells(2,1);
      pstring{1} = [pname1,' = ',parav1];
      pstring{2} = [pname2,' = ',parav2];
      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';
    if ~isempty(dud.plhand)    
      delete(dud.plhand);
    end    
    dud.plhand = [];    
   
    dud.color = sud.color;
    dud.fieldtype = Arrflag;
    dud.npts = NumbFPts;  

    set(dfdisp,'user',dud);
    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 = dwind - dud.settings.magn*[ud.DY;ud.DY];
    ud.cwind = wind - dud.settings.magn*[DY(1);-DY(1);DY(2);-DY(2)];
    ud.stop = 0;
    ud.gstop = 1;
    ud.plot = 1;
    set(dud.axes,'user',ud);
    ppkbd = findobj('name','DFIELD5 Keyboard input','vis','on');
    if ~isempty(ppkbd),dfield5('kbd'),end
    dfield5('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 DFIELD5 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 
      hhh = findobj('tag','zbmenu');
      set(hhh,'enable','on');
    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,'riceprint')
  
  printcode = [	
    'print             ';
    'print -Pbaker     ';
    'print -Pbrown     ';
    'print -Phanszen   ';
    'print -Pjones     ';
    'print -Plovett    ';
    'print -Psid       ';
    'print -Pwiess     ';
    'print -Pwrc       ';
    'print -Pfondren1ps';
    'print -Pfondren2ps';
    'print -Pmd105     '
    ];
  
  printname = [
    'The printer in this location';
    'Baker College               ';
    'Brown College               ';
    'Hanszen College             ';
    'Jones College               ';
    'Lovett College              ';
    'Sid Rich College            ';
    'Wiess College               ';
    'Will Rice College           ';
    'Fondren (front)             ';
    'Fondren (back)              ';
    'Mudd 105                    ';
    ];
  
  NN = size(printname,1);
  
  dud = get(gcf,'user');
  
  txtl = 150;txth = 21;
  
  pbwidth = 80;
  framewidth = txtl+pbwidth+20;
  frht = 10 + NN*txth +2*txth +10;
  frleft = 10;
  frbot = 20 + txth;
  frpos = [frleft frbot framewidth frht];
  figwidth = framewidth + 20;
  
  dfrice = figure('pos',[150, 50, figwidth, frbot+frht + 10],...
      'name','DFIELD5 Rice printer selection','numb','off',...
      'vis','off','tag','DFIELD5');
  
  dfield5('figdefault',dfrice);
  
  set(dfrice,'user',dud.printstr);		
  
  
  riceframe = uicontrol('style','frame','pos',frpos,'vis','off');
  printstr = dud.printstr;
  maxval = 1:NN;
  
  for j=1:NN
    
    rad(j) = uicontrol('style','radio',...
	'pos',[(frleft+8), frbot+10+(NN-j)*txth pbwidth+txtl txth],...
	'string',printname(j,:),'max',j,...
	'value',strcmp(printstr,deblank(printcode(j,:)))*maxval(j),...
	'tag',deblank(printcode(j,:)),'hor','left','vis','off');
  end    
  
  
  for i=1:NN
    set(rad(i),'user',rad(:,[1:(i-1),(i+1):NN]));
  end
  
  callrad = [
    'dfdisp = findobj(''name'',''DFIELD5 Display'');'...
	'dud = get(dfdisp,''user'');',...
	'me = get(gcf,''currentobject'');',...
	'kk = get(me,''max'');',...
	'set(get(me,''user''),''value'',0),',...
	'set(me,''value'',kk),',...
	'dud.printstr = get(me,''tag'');',...
	'set(dfdisp,''user'',dud);'];
  
  set(rad,'call',callrad);
  textpos = [frleft+8, frbot+10+(NN+0.5)*txth,pbwidth+txtl,txth];
  textstr = 'Select the printer of your choice';
  textw=uicontrol('style','text','pos',textpos,'string',textstr,...
      'hor','left','vis','off');
  bwidth=80;
  bpos = [(figwidth-bwidth)/2,10,bwidth,txth];
  offbutt = uicontrol('style','push','pos',bpos,'string','Go away',...
      'call','delete(gcf)');
  
  set(findobj(gcf,'type','uicontrol'),'units','normal');
  set(gcf,'resize','on');	
  set(gcf,'vis','on');
  set(get(gcf,'child'),'vis','on');

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','dfield5(''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')
    dfield5('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]);
  dfield5('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)];
    dfield5('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','dfield5(''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 DFIELD5 Keyboard figure.

  compute = 1;
  kh = get(gcf,'user');
  dfdisp = findobj('name','DFIELD5 Display');
  if (isempty(dfdisp))
    dfield5('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
    dfield5('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','dfield5(''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 DFIELD5 Keyboard figure which allows accurate input of 
  % initial values using the keyboard.
  
  ppkbd = findobj('name','DFIELD5 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','DFIELD5 Keyboard input',...
      'vis','off',...
      'numb','off','tag','DFIELD5');
  
  dfield5('figdefault',ppkbd);
  fmenu = findobj(gcf,'label','File');  % uimenu('label','File');
  clmenu = uimenu(fmenu,'label','Close',...
      'call','close','accelerator','W');
  compmenu = uimenu(fmenu,'label','Compute',...
      'call','dfield5(''kcompute'')','accelerator','C');
  qmenu = uimenu(fmenu,'label','Quit DFIELD5',...
      'call','dfield5(''quit'')','separator','on');
  
  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');
  
  kbd.t = uicontrol('style','text','string','<= t<= ');
  
  kbd.spec = uicontrol('style','check','horiz','center',...
      'string','Specify a computation interval.');
  
  kbd.comp = uicontrol('style','push',...		
      'string','Compute','call','dfield5(''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);
  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','dfield5(''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 = dwind - dud.settings.magn*[aud.DY;aud.DY];
  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','dfield5(''down'')',...
      'WindowButtonUpFcn','');
  dfield5('dirfield',disph);
  dfset = findobj('name','DFIELD5 Setup');
  if isempty(dfset)
    dfield5('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 solutions option.  
  
  disph = gcf;
  dud = get(disph,'user');
  notice = dud.notice;
  dud.eqpts = [];
  
  delete(dud.solhand);
  dud.solhand = [];
  delete(dud.plhand);
  dud.plhand = [];
  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','DFIELD5 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'',''DFIELD5 Display'');',...
	'dud = get(dfdisp,''user'');',...
	'dud.settings = data.settings;',...
	'set(dfdisp,''user'',dud);',...
	'dfset = findobj(''name'',''DFIELD5 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 = 1; nudge = 1; varl = 60;
	
  setfig = figure('name','DFIELD5 Solver settings','numb','off',...
      'tag','DFIELD5','vis','off');
  
  dfield5('figdefault',setfig);
  fmenu = findobj(gcf,'label','File');  % uimenu('label','File');
  clmenu = uimenu(fmenu,'label','Close',...
      'call','close',...
      'accelerator','W');

  qmenu = uimenu(fmenu,'label','Quit DFIELD5',...
      'call','dfield5(''quit'')',...
      'separator','on');
	
  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);'];
  
  
  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);		
  varl = varl*ht/19;
  setspfrb = 2*left + 2*ht;
  setspfrht = 3*nudge + 3*ht;
  setfrb = setspfrb + setspfrht;
  setfrh = 3*nudge + 3*ht;
  figh = setfrb + setfrh + left;
  gobb = left;
  chbb = gobb + ht;
  rtolb = setfrb + nudge;
  d1b = rtolb + ht;
  ssizeb = (rtolb+d1b)/2;
  ssb = d1b + ht + nudge;
  ssw = max(max(ext1(3) + varl,ext2(3)),200);
  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,figw,setspfrht]);
  set(data.speed,'pos',[sl,spb1,ssw-sl,ht]);
  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,rtolb,stw,ht],...
      'string','Relative error tolerance: '); 
      
    data.rtol = uicontrol('style','edit','call',rtolcall,...
      'pos',[sl+stw,rtolb,varl,ht],...
      'string',num2str(data.settings.tol));
  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 DFIELD5 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','DFIELD5 Window settings');
  if ~isempty(setfig)
    close(setfig);
  end
  setfig = figure('name','DFIELD5 Window settings','numb','off',...
      'tag','DFIELD5','vis','off');
  
  dfield5('figdefault',setfig);
  fmenu = findobj(gcf,'label','File');  % uimenu('label','File');
  clmenu = uimenu(fmenu,'label','Close',...
      'call','close',...
      'accelerator','W');
  chmenu = uimenu(fmenu,'label','Change Settings',...
      'call','dfield5(''setchange'')',...
      'accelerator','C');
  qmenu = uimenu(fmenu,'label','Quit DFIELD5',...
      'call','dfield5(''quit'')',...
      'separator','on');
	

  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','dfield5(''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 DFIELD5 Settings window.

  dfsett = gcf;
  data = get(dfsett,'user');
  dfdisp = findobj('name','DFIELD5 Display');
  dfset = findobj('name','DFIELD5 Setup');
  if isempty(dfdisp)|isempty(dfset)
    dfield5('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);
    dfield5('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,'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,'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];
  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)';     
    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','dfield5(''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 DFIELD5 Display Window.  Please be patient.';
  set(dud.notice,'string',nstr);
  
  nstr(1:4) = nstr(2:5);
  nstr{5} = 'Printing the DFIELD5 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 = cells(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',[400,200],...
      'initialvalue',wch,...
      'name','DFIELD5 Zoomback',...
      'promptstring','Select a rectangle:',...
      'OKString','Zoom');
  
  if (~isempty(sel))
    WINvect = wmat(sel,:);
    dud.syst.wind = WINvect;
    set(gcf,'user',dud);
    dfset = findobj('name','DFIELD5 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);
    dfield5('dirfield',disph);
  end
  
elseif strcmp(action,'figdefault')

  fig = input1;
  set(fig,'CloseRequestFcn','dfield5(''closefcn'')');
  dfset = findobj('name','DFIELD5 Setup');
  sud = get(dfset,'user');
  fs = sud.fontsize;
  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 'display'
    whitebg(fig,[1 1 1])
    set(fig,'defaultaxescolor',[1 1 1])
    set(fig,'defaultaxescolororder', ...
        [0 0 1;0 .5 0;1 0 0;0 .75 .75;.75 0 .75;.75 .75 0;.25 .25 .25]) % bgrymck
    set(fig,'colormap',jet(64))
    set(fig,'defaultsurfaceedgecolor',[0 0 0])
    set(fig,'color',[1 1 1]*240/255);
    set(fig,'defaultuicontrolbackgroundcolor',[1 1 1]*220/255);

  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,'menubar','none');
  uimenu(fig,'label','File',...
      'Handle','on');
  name = get(fig,'name');
  if strcmp(name,'DFIELD5 Setup') | strcmp(name,'DFIELD5 Display')...
	| strcmp(name,'DFIELD5 Linearization')
    uimenu(fig,'label','Edit',...
	'Handle','on');
  end
  wind = uimenu(fig,'label','Window',...
      'Handle','on');
  hhsetup = get(0,'ShowHiddenHandles');
  set(0,'ShowHiddenHandles','on');
  wh = findobj('type','uimenu','label','Window');
  fh = get(0,'children');
  fh(find(fh == fig)) = [];
  for j=1:length(wh)
    uimenu(wh(j),'label',name,...
	'accelerator',int2str(fig),...	
	'user',fig,...
	'call','figure(get(gcbo,''user''))');
  end
  for j=1:length(fh)
    name = get(fh(j),'name');
    if strcmp(name,'')
      name = ['Figure No. ',int2str(fh(j))];
    end  
    uimenu(wind,'label',name,...
	'accelerator',int2str(fh(j)),...		
	'user',fh(j),...
	'call','figure(get(gcbo,''user''))');
  end
  set(0,'ShowHiddenHandles',hhsetup);
  
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')
  
  dfdisp = findobj('name','DFIELD5 Display');
  if (~isempty(dfdisp))
    dud = get(dfdisp,'user');
    fcn = dud.function;
    if (exist(fcn)==2) delete([fcn,'.m']);end
  end
  h = findobj('tag','DFIELD5');
  delete(h);
  
elseif strcmp(action,'restart')
  
  dfdisp = findobj('name','DFIELD5 Display');
  dfset = findobj('name','DFIELD5 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','DFIELD5');
  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,'DFIELD5 Setup') | strcmp(name,'DFIELD5 Display')
    quest = ['Closing this window will cause all DFIELD5 ',...
	  'windows to close, and DFIELD5 will stop.  ',...
	  'Do you want to quit DFIELD5?'];
    butt = questdlg(quest,'Quit DFIELD5?','Quit','Cancel','Quit');
    if strcmp(butt,'Quit')
      dfield5('quit');
    end
  elseif strcmp(name,'DFIELD5 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,'confused')
  
  tstring = 'DFIELD5 is totally confused';
  qstring = {['You will have to restart DFIELD5 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 DFIELD5.';
  bstr2 = 'Just quit DFIELD5.';
  bstr3 = 'Do nothing.';
  answer = questdlg(qstring,tstring,bstr1,bstr2,bstr3,bstr1);
  if strcmp(answer,bstr1)
    delete(findobj('tag','DFIELD5'));
    dfield5;return
  elseif strcmp(answer,bstr2)
    delete(findobj('tag','DFIELD5'));
    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;
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',30,'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',30,'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',30,'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',30,'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








