%
% Check Van Der Pol equation
% this is a stiff test
%

% set some option for the ODE numerical solver
options = odeset('RelTol',1e-4) ; %,'AbsTol',[1e-7 1e-7 1e-7]);

% se initial condition for the problem
Y0 = [ 2 0 ] ;

% Set a vector specifying the interval of integration, [t0,tf].
% The solver imposes the initial conditions at TSPAN(1),
% and integrates from TSPAN(1) to TSPAN(end). 
% To obtain solutions at specific times (all increasing or all decreasing), 
% use TSPAN = [t0,t1,...,tf].
TSPAN = [ 0:1:3000 ] ;

% load the "pointer" of the function to FUN
FUN = @vanderpol ;

% This is a stiff problem, the first three solver do not work! 
%[T1,Y1] = ode23   ( FUN, TSPAN, Y0 ) ;
%[T2,Y2] = ode45   ( FUN, TSPAN, Y0 ) ;
%[T3,Y3] = ode113  ( FUN, TSPAN, Y0 ) ;

% those solvers can solve a stiff problem
[T4,Y4] = ode15s  ( FUN, TSPAN, Y0 ) ;
[T5,Y5] = ode23s  ( FUN, TSPAN, Y0 ) ;
[T6,Y6] = ode23t  ( FUN, TSPAN, Y0 ) ;
[T7,Y7] = ode23tb ( FUN, TSPAN, Y0 ) ;

%
% plot the functions approximated by various solver
%
% first component
plot( T4, Y4(:,1), 'c-.', ...
      T5, Y5(:,1), 'm--',  ...
      T6, Y6(:,1), 'y:',  ...
      T7, Y7(:,1), 'k-.' )
      
% create a new window to plot
figure()

% second component
plot( T4, Y4(:,2), 'c-.', ...
      T5, Y5(:,2), 'm--',  ...
      T6, Y6(:,2), 'y:',  ...
      T7, Y7(:,2), 'k-.' )

