dfield7.m

matlab of linear system

	else
	  eulch = 'off';
	end
	rk2set.refine = 1;
	rk2set.tol = 0;
	rk2set.solver = 'Runge-Kutta 2';
	rk2set.solvhandle = @dfrk2;
	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.solvhandle = @dfrk4;
	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','dfield7(''settings'')');
	
	menu(6) = uimenu(menu(1),'label','Make the Display Window inactive.',...
			 'call','dfield7(''hotcold'')','separator','on');
	
	% View menu
	
	set(findobj(gcf,'label','&View'),'pos',4);
	set(findobj(gcf,'label','&Figure Toolbar'),...
	    'call','dfield7(''showbar'')');
	dud.menu = menu;
	
	set(dfdisp,'ToolBar','none');
	set(0,'showhiddenhandles',hhsetup);
	set(gcf,'WindowButtonDownFcn','dfield7(''down'')');
	set(dfdisp,'WindowButtonMotionFcn','dfield7(''cdisp'')');
	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);
      pf = pflag;
      k = find(pflag);
      if ~isempty(k)
	tstr = [tstr,'             '];
	for j = 1:length(k)
	  tstr = [tstr,'    ',pstring{k(j)}];
	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','dfield7 Keyboard input','vis','on');
	if ~isempty(ppkbd),dfield7('kbd'),end
	dfield7('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 dfield7 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','dfield7(''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')
      dfield7('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]);
    dfield7('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]=dfginput;
    NN = length(X);
    for k = 1:NN
      initpt = [X(k),Y(k)];
      dfield7('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','dfield7(''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');
      
      solvhandle = settings.solvhandle;
      
      cflag = 0;
      if intplus(2)>intplus(1)
	cflag = cflag + 1;
	if notice
	  nstr = get(notice,'string');