function [ ts, ys ] = matlabFailure(show_failure, tmax)
% matlabFailure: Demonstrate an instance of ode45 failure.
%
%   [ ts, ys ] = matlabFailure(tmax)
%
% Demonstrates the coupled oscillator example from
%
%   "Analysis Still Matters: A Surprising Instance of Failure of
%    Runge-Kutta-Felberg ODE Solvers", Joseph Skufca,
%    SIAM Review, v.46, n.4, pp.729-737.
%
% The oscillators should synchronize and stay synchronized, but 
%   in practice they synchronize and then destabilize around t = 100.
%
% The failure occurs because the true solution is linear, and so the
%   error in order p, p+1 both go to zero.  The integrator responds
%   by increasing the stepsize until instability sets in.
%
% In this particular case, the growth of the solution allows the
%   error tolerance to grow larger, and eventually the oscillator
%   pair slips into synchronization at another offset.  The effect
%   continues to grow since the error bound grows.
%
% Input Parameters:
%
%   show_failure:  Boolean.  Show the failure or not?  Default = 1. 
%
%   tmax:  Double.  Final time of integration.  Default = 200.
%
% Output Parameters:
%
%   ts: Vector.  Timesteps chosen by the integrator (column vector).
%
%   ys: Matrix.  Solution at each timestep, as returned by the integrator.
%
% Ian Mitchell for CS 403, 2/11/05.

% Check the input arguments.
if(nargin < 1)
  show_failure = 1;
end

if(nargin < 2)
  tmax = 200;
end

% Problem parameters.
omegas = [ 1; 1.5 ];
ks = [ 1; 1 ];
yinit = [ 3; 0 ];
options = odeset('Refine', 1);

% Simulate.
if(show_failure)
  % The standard ode45 desynchronizes after t = 100.
  [ ts, ys ] = ode45(@oscillators, [ 0, tmax ], yinit, options, omegas, ks);
else
  % The stiff solver happens to work, although not because the system is
  % stiff.  For more details, see 
  [ ts, ys ] = ode15s(@oscillators, [ 0, tmax ], yinit, options, omegas, ks);
end

% Plot the trajectories.
figure;
plot(ts, ys);

% Plot the timesteps.
figure;
plot(ts(1:end-1), diff(ts));



%---------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------------------------------------------------------------
function ydot = oscillators(t, y, omegas, ks)
% oscillators: implements the coupled oscillator ODE.
%
%   ydot =  oscillators(t, y, omegas, ks)
%
% Implements the coupled oscillator system from Skufca.
%
%    \dot y_1 = \omega_1 + k_1 \sin(y_2 - y_1)
%    \dot y_2 = \omega_2 + k_2 \sin(y_1 - y_2)
%
% Parameters:
%
%   t        Current time.
%   y        Current state.
%   omegas   Natural frequencies of the two oscillators (column vector).
%   ks       Coupling constants of the two oscillators (column vector).
%
%   ydot     State derivative.

delta = y(2) - y(1);

ydot = omegas + ks .* sin([ delta; -delta ]);