function vdpeig(action)

% This demo shows some interesting details about the solution of the VanDer Pol 
% second order ODE. The solution is drawn in the top of the picture with a dot 
% moving around dragging a tail of the last 40 steps. 
% At the bottom is drawn a picture of the eigenvalues in the complex plane, 
% the scale is {-10 .. 10] in both directions, the origin is indicated by a green star.
% Together with the eigenvalues the points scaled by the stepsize is drawn as
% a green dot, it moves near the negative real axis and the endpoint of the stability region 
% is indicated by a green cross.
% When the demo runs the points will move around their respective closed curves illustrating
% what happens in the stepsize control, how it uses the stability interval and how the tail 
% gets short when the stiffness ratio gets large.
% The GUI is just the one adapted for the ODEdemo's in the Matlab distribution.
% The demo uses the two files vdpol.m and vdpjac.m
% The demo hs been tested in versions 5.3.1 (Unix) and student edition 4a (Mac)
% Submitted by Per Grove Thomsen
% Institute for Mathematical Modelling
% Danish Technical University
% Building 305, DK 2800 Lyngby , Denmark.
% E-mail: pgt@imm.dtu.dk
% URL: http://www.imm.dtu.dk/documents/users/pgt/homepage.html

play= 1;
stop=-1;

if nargin<1,
    action='initialize';
end;

if strcmp(action,'initialize'),
    oldFigNumber=watchon;

    figNumber=figure( ...
        'Name','Van Der Pol oscillator', ...
        'NumberTitle','off', ...
	'Visible','off');
    axes( ...
        'Units','normalized', ...
        'Position',[0.05 0.10 0.75 0.95], ...
        'Visible','off');

    text(0,0,'Press the "Start" button to see the VanDerPol demo', ...
	'HorizontalAlignment','center');
    axis([-1 1 -1 1]);

    %===================================
    % Information for all buttons
    labelColor=[0.8 0.8 0.8];
    yInitPos=0.90;
    xPos=0.85;
    btnLen=0.10;
    btnWid=0.10;
    % Spacing between the button and the next command's label
    spacing=0.05;

    %====================================
    % The CONSOLE frame
    frmBorder=0.02;
    yPos=0.05-frmBorder;
    frmPos=[xPos-frmBorder yPos btnLen+2*frmBorder 0.9+2*frmBorder];
    h=uicontrol( ...
        'Style','frame', ...
        'Units','normalized', ...
        'Position',frmPos, ...
	'BackgroundColor',[0.50 0.50 0.50]);
    
    %====================================
    % The START button
    btnNumber=1;
    yPos=0.90-(btnNumber-1)*(btnWid+spacing);
    labelStr='Start';
    cmdStr='start';
    callbackStr='vdpeig(''start'');';

    % Generic button information
    btnPos=[xPos yPos-spacing btnLen btnWid];
    startHndl=uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'Interruptible','on', ...
        'Callback',callbackStr);
   
    %====================================
    % The STOP button
    btnNumber=3;
    yPos=0.90-(btnNumber-1)*(btnWid+spacing);
    labelStr='Stop';
    % Setting userdata to -1 (=stop) will stop the demo.
    callbackStr='set(gca,''Userdata'',-1)';
    
    % Generic  button information
    btnPos=[xPos yPos-spacing btnLen btnWid];
    stopHndl=uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'Enable','off', ...
        'String',labelStr, ...
        'Callback',callbackStr);

    %====================================
    % The INFO button
    labelStr='Info';
    callbackStr='vdpeig(''info'')';
    infoHndl=uicontrol( ...
        'Style','push', ...
        'Units','normalized', ...
        'position',[xPos 0.20 btnLen 0.10], ...
        'string',labelStr, ...
        'call',callbackStr);

    %====================================
    % The CLOSE button
    labelStr='Close';
    callbackStr='close(gcf)';
    closeHndl=uicontrol( ...
        'Style','push', ...
        'Units','normalized', ...
        'position',[xPos 0.05 btnLen 0.10], ...
        'string',labelStr, ...
        'call',callbackStr);
    
    % Uncover the figure
    hndlList=[startHndl stopHndl infoHndl closeHndl];
    set(figNumber,'Visible','on', ...
	'UserData',hndlList);

    watchoff(oldFigNumber);
    figure(figNumber);

elseif strcmp(action,'start'),
    axHndl=gca;
    figNumber=gcf;
    hndlList=get(figNumber,'UserData');
    startHndl=hndlList(1);
    stopHndl=hndlList(2);
    infoHndl=hndlList(3);
    closeHndl=hndlList(4);
    set([startHndl closeHndl infoHndl],'Enable','off');
    set(stopHndl,'Enable','on');

    % ====== Start of Demo
    set(figNumber,'Backingstore','off');
    % The graphics axis limits are set to values known 
    % to contain the solution.
    set(axHndl, ...
        'XLim',[-10 10],'YLim',[-10 10],'ZLim',[0 10], ...
	'Userdata',play, ...
        'XTick',[0],'YTick',[0],'ZTick',[7], ...
	'Drawmode','fast', ...
        'Visible','on', ...
	'NextPlot','add', ...
        'Userdata',play, ...
	'View',[-37.5,30]);
    xlabel('y1');
    ylabel('y2');
    zlabel('Z');
    vdpstr1 = ['Van Der Pol problem                 '];
    vdpstr2 = ['y" = -mu (1-y^2)y` - y              '];
    vdpstr3 = ['for  mu = 5.0                       '];
    %text(-9,-9,5,vdpstr,'fontsize',12);

    % The values of the global parameters are
    global mu
    mu = 5.;
   
    % The numerical integration, and the display of the evolving solution,
    % are handled by the function ODE23P.
 
    FunFcn='vdpol';
    % The initial conditions below will produce good results
    %y0 = [ 5 0];
    
    y0(1)=2;
    y0(2)=0;
   
    t0=0;
    tfinal=100;
    pow = 1/3;
    tol = 0.0001;
 
    t = t0;
    hmax = (tfinal - t)/5;
    hmin = (tfinal - t)/20000;
    h = (tfinal - t)/100;
    y = y0(:);
    yd=[1.5*y; 7];
    tau = tol * max(norm(y,'inf'),1);
    % Compute the eigenvalues of the Jacobian.
    z=vdpjac(t,y,mu);
            e=eig(z);
            al1=real(e(1)); be1 = imag(e(1));
            al2=real(e(2)); be2 = imag(e(2));
            ev1 =[al1 be1]';
            ev2 =[al2 be2]';
    % Save L steps and plot like a comet tail.
    L = 40;
    Y = yd*ones(1,L);
  	E1 = ev1*ones(1,L);
 	E2 = ev2*ones(1,L);
 	hE1= h*E1;
 	hE2= h*E2;
  	zz = [0;0;0];
  	stb =[-2.5; 0;0];
    cla;
   % figure(1)
   nulpkt = line( ...
	'color','g', ...
	'linestyle','*', ...
	'markersize',15, ...
	'erase','none', ...
	'xdata',zz(1),'ydata',zz(2),'zdata',zz(3));
	nulpkt = line( ...
	'color','g', ...
	'linestyle','+', ...
	'markersize',15, ...
	'erase','none', ...
	'xdata',stb(1),'ydata',stb(2),'zdata',stb(3));
    head = line( ...
	'color','r', ...
	'linestyle','.', ...
	'markersize',25, ...
	'erase','xor', ...
	'xdata',yd(1),'ydata',yd(2),'zdata',yd(3));
    body = line( ...
	'color','y', ...
	'linestyle','-', ...
	'markersize',5, ...
	'erase','none', ...
	'xdata',[],'ydata',[],'zdata',[]);
    tail=line( ...
	'color','g', ...
	'linestyle','-', ...
	'markersize',5, ...
	'erase','none', ...
	'xdata',[],'ydata',[],'zdata',[]);
	scaled1 = line( ...
	'color','g', ...
	'linestyle','.', ...
	'markersize',15, ...
	'erase','none', ...
	'xdata',hE1(1),'ydata',hE1(2));
	scaled2 = line( ...
	'color','g', ...
	'linestyle','.', ...
	'markersize',15, ...
	'erase','none', ...
	'xdata',hE2(1),'ydata',hE2(2));
	
	text(-19,-9,10,vdpstr1,'fontsize',12);
	text(-19,-9,9,vdpstr2,'fontsize',12);
	text(-19,-9,8,vdpstr3,'fontsize',12);
	
% and now the eigenvalues
%	figure(2)
	head1 = line( ...
	'color','r', ...
	'linestyle','+', ...
	'markersize',10, ...
	'erase','xor', ...
	'xdata',ev1(1),'ydata',ev1(2));
    body1 = line( ...
	'color','r', ...
	'linestyle','.', ...
	'markersize',15, ...
	'erase','none', ...
	'xdata',[],'ydata',[]);	
	tail1 = line( ...
	'color','r', ...
	'linestyle','o', ...
	'markersize',5, ...
	'erase','none', ...
	'xdata',[],'ydata',[]);
	head2 = line( ...
	'color','b', ...
	'linestyle','o', ...
	'markersize',10, ...
	'erase','xor', ...
	'xdata',ev2(1),'ydata',ev2(2));
    body2 = line( ...
	'color','b', ...
	'linestyle','.', ...
	'markersize',10, ...
	'erase','none', ...
	'xdata',[],'ydata',[]);	
	tail2 = line( ...
	'color','g', ...
	'linestyle','o', ...
	'markersize',5, ...
	'erase','none', ...
	'xdata',[],'ydata',[]);
 
    % The main loop
    while (get(axHndl,'Userdata')==play) & (h >= hmin)
   	if t + h > tfinal, h = tfinal - t; end
   	% Compute the slopes
      	s1 = feval(FunFcn, t, y,mu);
      	s2 = feval(FunFcn, t+h, y+h*s1,mu);
      	s3 = feval(FunFcn, t+h/2, y+h*(s1+s2)/4,mu);
 
      	% Estimate the error and the acceptable error
      	delta = norm(h*(s1 - 2*s3 + s2)/3,'inf');
      	tau = tol*max(norm(y,'inf'),1.0);
 
      	% Update the solution only if the error is acceptable
      	ts = t;
      	ys = y;
      	if delta <= tau
            t = t + h;
            y = y + h*(s1 + 4*s3 + s2)/6;
            yd=[1.5*y ;7];
         % Compute the eigenvalues of the Jacobian.   
            z=vdpjac(t,y,mu);
            e=eig(z);
            al1=real(e(1)); be1 = imag(e(1));
            al2=real(e(2)); be2 = imag(e(2));
 		   	ev1=[al1,be1]';
 			   ev2=[al2,be2]';
 		   	hE1 = h*ev1;
 		   	hE2 = h*ev2;
            % Update the plot
            Y = [yd Y(:,1:L-1)];
            E1 = [ev1 E1(:,1:L-1)];
            E2 =[ev2 E2(:,1:L-1)];
            set(head,'xdata',Y(1,1),'ydata',Y(2,1),'zdata',Y(3,1))
            set(body,'xdata',Y(1,1:2),'ydata',Y(2,1:2),'zdata',Y(3,1:2))
            set(tail,'xdata',Y(1,L-1:L),'ydata',Y(2,L-1:L),'zdata',Y(3,L-1:L))
            drawnow;
            set(head2,'xdata',E2(1,1),'ydata',E2(2,1))
            set(body2,'xdata',E2(1,1:2),'ydata',E2(2,1:2))
            drawnow;
    		   set(head1,'xdata',E1(1,1),'ydata',E1(2,1))
            set(body1,'xdata',E1(1,1:2),'ydata',E1(2,1:2))
            set(scaled1,'xdata',hE1(1,1),'ydata',hE1(2,1))
            set(scaled2,'xdata',hE2(1,1),'ydata',hE2(2,1))
            drawnow;
        end
 
      	% Update the step size
      	if delta ~= 0.0
            h = min(hmax, 0.9*h*(tau/delta)^pow);
      	end

    end;    % Main loop ...
    % ====== End of Demo
    set([startHndl closeHndl infoHndl],'Enable','on');
    set(stopHndl,'Enable','off');

elseif strcmp(action,'info');
    ttlStr='Van Der Pol Info';
    hlpStr= ...                                             
        ['                                               '  
         ' This demo animates the integration of  the    '  
         ' two coupled nonlinear differential equations  '  
         ' that define the "Van Der Pol oscillator", a   '  
         ' system first described by Duch engineer B.    '
         ' Van Der Pol in the 1920`s                     '  
         ' A model of an electrical circuit that is used '
         ' to generate stable oscillating signals.       '  
         '                                               '  
         ' As the integration proceeds you will see a    '  
         ' point moving around in a stable orbit in      '  
         ' 2-D space . The eigenvalues of the Jacobian   '
         ' are computed and displaued in the complex     '
         ' plane. They also move around in a closed graph'
         ' Mu is equal to 5 .                            '
         'The orbit is bounded,periodic and convergent . '  
         '                                               '  
         ' Use the "Start" and "Stop" buttons to control '  
         ' the animation.                                '  
         '                                               '  
         ' File names: vdpol.m  vdpeig.m  vdpjac.m       '];
      
    helpfun(ttlStr,hlpStr);                                 

end;    % if strcmp(action, ...

