% This code uses matlab's ode45 to numerically solve the van der
% Pol equation and uses the event finder to locate zeros, thereby
% estimating the period.
% This process is performed for various values of $\mu$ and the
% results dumped into a file.
% Use vdp_period_plot to load that file and process the data.

% 2003-11-14 cbm
% no rights reserved: public domain code


% where to dump the results
filename = 'vdp_period_datafile2';

% the mu values to compute over
muset = logspace(log10(1), log10(100), 100);

% integrate from 0 to tf for each mu
tf = 500;

% initial condition for each integration
x0 = [1; -2/3];



% matlab can make you lazy about initializing arrays.  If you
% choose not to then codes can be difficult to debug (unless you
% put a "clear" at the top of your file
periods = zeros(size(muset));

% see 'help ode45' for more information on events
opt = odeset('Events', @vdp_period_events);
% defaults were not accurate enough (perhaps events are found to a
% lesser accuracy?)
opt = odeset(opt, 'RelTol', 1e-5, 'AbsTol', 1e-8);

for m=1:length(muset)
  mu = muset(m)

  [t, x, te, xe] = ode45(@vdp_period_f, [0 tf], x0, opt, mu);
  % te: times when events (roots in this case) occured
  % xe: values at each events (xe(1,:) ~= 0)

  
  disp(['detected '  num2str(length(te))  ' events']);

  % throw away this many roots (maybe we don't trust the first one
  % in case it occurs during the transient)
  TA = 1;

  % time between events
  tbe = te((2+TA):end) - te((1+TA):(end-1));

  % there are two events / period
  periods(m) = 2*mean(tbe);
  
  disp(['period estimate: '  num2str(periods(m))  '']);

  figure(1);  clf;  hold on;
  axis([-6 6 -10 10]);
  plot(x0(1), x0(2), 'ko');
  plot(x(:,1), x(:,2), 'bx');
  H = plot(x(:,1), x(:,2), 'k-');
  set(H, 'color', .85*[1 1 1]);
  xlabel('x(t)-axis');
  ylabel('x''(t)-axis');   % note '' will render as '
  title('phase portrait');

  figure(2);  clf;  hold on;
  plot([0 tf], [0 0], 'k:');
  plot(t, x(:,1), 'bx-');
  plot(te, xe(:,1), 'r*');
  title(['van der Pol solution: \mu = '  num2str(mu)  ', period ~= ' ...
	  num2str(periods(m))  ]);
  xlabel('t-axis');
  ylabel('x-axis');
  legend('x = 0', 'numerical solutions', 'events (roots)', 4);
  print('-dpng', '-r40', ['vdp_solution_mu'  num2str(mu)  '.png']);
  
  pause(0);
end

save(filename, 'muset', 'periods');


