% Consider the second order differential equation known as the
% van der Pol equation
%             y'' + mu*(y^2 - 1)y' + omegasq*y = 0
% where the prime ' denotes the time-derivative operator.  We can rewrite
% this as a system of first order differential equations:
%             y1' = y2
%             y2' = mu*(1 - y1^2)-omegasq*y2
%
% To simulate the above system, we create an m-file vdpsim.m that calls
% a function M-file vdpol.m which returns the state derivative vector. It invokes ode23.m
% vdpsim.m plots time responses of y and ydot in Figure 1 and the phase plane plot in Figure 2  
% It also queries new inputs for the parameter mu and the initial condition y(0)

clear all;

status='y';
mu=1.0; omegasq=10.;            % default values of parameters
tspan=[0.0 20.0];               % initial time t0 = 0; final time tfinal = 20;
y0 = [0.1 0.0]';                % default initial conditions y(0)=0  ydot(0)=0.1;

disp('    ')
disp('Van der Pol Equation  d2y/dt2 + mu*(y^2 - 1)dy/dt + omegasq*y = 0')
disp('    ')

while (status=='y')
   mu=input('enter the value of mu (default mu=1.0) ');
   y0(1)=input('enter the value of y0(1) (default y0(1)=0.1) ');
   options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]);  %options=[];
   [t,y] = ode23('vdpol',tspan,y0,options,mu,omegasq);

   figure(1)
   plot(t,y) 
   title('van der Pol equation time history')
   xlabel('time t')
   ylabel('position y and velocity ydot')
   grid on
   
   figure(2)
   plot(y(:,1),y(:,2)) 
   title('van der Pol equation - phase plane plot')
   xlabel('position y')
   ylabel('velocity ydot')
   grid on
   
   status=input('enter y to continue or anything else to terminate  ','s')
end