clear all; close all;

% Parameters
tspan = [0 3000];
y0 = [2 0];
mu = 1000;

% right hand side of the Van der Pol ODE
% (already written as a first order system :-)
f = @(t,y) [y(2); mu*(1. - y(1)^2)*y(2) - y(1)];

% Jacobi matrix of the right hand side of Van der Pol ODE
J = @(t,y) [0,1;-2.*mu*y(1)*y(2) - 1,mu*(1. - y(1)^2);];

% explicit solver: ode45
disp('Solver: ode45');
tic;
[t45 ,y45 ] = ode45 (f,tspan,y0);
toc
% compute stiffness ratio S along the solution curve, i.e.
% for each time step
S45  = zeros(size(t45 ));
for istep=1:length(t45 )
  % compute eigenvalues of Jacobi matrix (simply with eig)
  lmbdas = eig(J(t45 (istep),y45 (istep,:)));
  % compute stiffness ratio
  S45 (istep) = max(abs(lmbdas))/min(abs(lmbdas));
end

% implicit solver: ode23s
disp('Solver: ode23s');
tic;
[t23s,y23s] = ode23s(f,tspan,y0);
toc
% compute stiffness ratio S along the solution curve, i.e.
% for each time step
S23s = zeros(size(t23s));
for istep=1:length(t23s)
  % compute eigenvalues of Jacobi matrix (simply with eig)
  lmbdas = eig(J(t23s(istep),y23s(istep,:)));
  % compute stiffness ratio
  S23s(istep) = max(abs(lmbdas))/min(abs(lmbdas));
end

% plot solution
figure(1);
plot(t45 ,y45 (:,1),'r--x')
hold on;
plot(t23s,y23s(:,1),'b-o')
title('Van der Pol ODE (solution y)')
xlabel('t');
ylabel('y')
legend('ode45 (explicit)','ode23s (implicit)');

% plot step size for each step
figure(2);
semilogy(t45 (2:end),diff(t45 ),'b-x');
hold on;
semilogy(t23s(2:end),diff(t23s),'r-o');
title('Van der Pol ODE (step size h)')
xlabel('t');
ylabel('h')
legend('ode45 (explicit)','ode23s (implicit)');

% plot stiffness ratio along solution curve
figure(3)
semilogy(t23s,S23s,'b-o')
title('Van der Pol ODE (stiffness ratio S)')
xlabel('t');
ylabel('S')

