function vanderpol(alpha,y,z,Fmx)

%  Phase portrait and Poincare map for Van der Pol problem

% Inputs; 
%  y, z   vectors of values defining all coordinate pairs v=(y,z)
%         at which to plot the vector F = (dy/dt,dz/dt). The 
%         positive part of the plotted range of y is also used to make
%         a Poincare map
%  alpha  damping parameter
%  s
% Van der Pol derivative function is at end of script.  We use
% ode15s rather than ode45 in this script because if alpha>>1 the relaxation
% oscillation can involve severe changes in the required adaptive
% timestep to follow the entire orbit accurately.
  
% Plotting parameters (s for scaling vector lengths, Fmx for
% setting a max vector length).

  s = 0.4;
  
  subplot(2,2,1)
  PhasePortrait(@dudt,y,z,s,Fmx)
  xlabel('y')
  ylabel('z = dy/dt')
  title(['Van der Pol oscillator phase portrait for \alpha = ' num2str(alpha)])

  % Plot unstable fixed point as open black circle
  
  hold on
  plot(0,0,'ko')
  
  % Plot an orbit in red. Orbit parameters are plotting frequency
  % dt, initial position u0, and max time shown tmax.

  dt = 0.1;
  u0 = [0.5; 0.5];
  tmax = max(3*alpha,10);

  tspan = 0:dt:tmax;
  [t, u] = ode15s(@dudt,tspan,u0);
  plot(u(:,1),u(:,2),'r.')
  hold off  

  subplot(2,2,2)
  
  % Poincare map along z=0, y>0 using event detection option of ode45
  
  options = odeset('Events',@events);

  n0 = 25;
  max(y)
  y0 = max(y)*(1:n0)/n0;   % Initial y values (z0=0)
  y1 = zeros(1,n0); % Vector to be filled with corresponding y after one orbit 
  for j = 1:n0
    tspan = [0 tmax];
    u0 = [y0(j); 0];
    [t,u,tev,ue,ie] = ode15s(@dudt,tspan,u0,options);
    y1(j) = ue(end,1);
  end
  plot(y0,y1,'.',y0,y0,'--')
  xlabel('y_n')
  ylabel('y_{n+1}')
  legend('y_{n+1}','1:1 line','Location','SouthEast')
  title(['Van der Pol Poincare map for \alpha = ' num2str(alpha)])
  
%-----------------------------------------------------------------  
function dyzdt = dudt(t,u)
  dyzdt = [u(2,:); alpha*(1-u(1,:).^2).*u(2,:) - u(1,:)]; 
end   % End nested function dudt
%-----------------------------------------------------------------
function [value,isterminal,direction] = events(t,u)

  value = u(2);  % Detect if z=0
  isterminal = 1; % Stop integration 
  direction = -1; % ...but only if z is decreasing (to ensure y>0).
end   % End nested function events

end