% runvdp(mu)
% Demonstrate relative performance of ode45 and ode15s for
% the van der Pol oscillator as a function of mu.
%
function runvdp(mu)

tau = (3-2*log(2))*mu + 4.676*mu^(-1/3); % Estimated period
tspan = [0, 3*tau]; % Go about 3 cycles
y0 = [2; 0]; % Initial conditions
opt = odeset('Stats', 'on'); % Print diagnostics
ode = @(t,y) vanderpoldemo(t,y,mu); % MATLAB already has f_vdp

fprintf('\n--- ODE45 solve ---\n');
tic; [tn,yn] = ode45(ode, tspan, y0, opt); toc

fprintf('\n--- ODE15s solve ---\n');
tic; [ts,ys] = ode15s(ode, tspan, y0, opt); toc

figure(1);
plot(tn,yn(:,1), 'r-', ts,ys(:,1), 'b-');
xlabel('t');
ylabel('x');
title(sprintf('Van der Pol time history, mu = %f', mu));

figure(2);
plot(yn(:,1), yn(:,2), 'r-', ys(:,1), ys(:,2), 'b-');
xlabel('x');
ylabel('y');
title(sprintf('Van der Pol phase plane, mu = %f', mu));