function duffing ( tspan, delta, iplot )
% Duffing
%
% This m-file is part of the material of the lectures on
%   Dynamical Systems and Nonlinear Analysis
% at the Math Department of Halle University, Germany.
%
% Author :      Martin Arnold, martin.arnold@mathematik.uni-halle.de
% Version of :  Oct  8, 2009
%
% Time integration of the Duffing equation
%
% Parameters:
%   tspan  (input)  : time span
%   delta  (input)  : parameter in Duffing equation
%   iplot  (input)  : control flag output
%                       =1 .. (t,u) plot
%                       =2 .. phase plot
%
% Example:
%   duffing ( [ 0,  10 ], 0.1,  1 );
%   duffing ( [ 0, 100 ], 0.01, 2 );

% -> define initial values, rhs etc.
a = 0.0;
b = 1.0;

rhs = @(t,u) ( duffrhs ( t, u, delta ) );

% -> time integration
[ t, u ] = ode45 ( rhs, tspan, [ a b ]' );

% -> visualization
if (iplot==1),

    plot ( t, u );

    legend ( 'x(t)', 'x''(t)' );

    xlabel ( 't' );
    ylabel ( 'u' );

elseif (iplot==2),

    plot ( u(:,1), u(:,2) );
    hold on;
    plot ( u(1,1), u(1,2), '.', 'MarkerSize', 24 );
    hold off;

    xlabel ( 'x' );
    ylabel ( 'x''' );

    set ( gca, 'DataAspectRatio',    [ 1 1 1 ] );
    set ( gca, 'PlotBoxAspectRatio', [ 1 1 1 ] );

end;

title ( [ 'Duffing equation.  Parameter  \delta = ', sprintf('%0.4g',delta) ] );

% ----------------------------------------------------------------- %

function uprime = duffrhs ( t, u, delta )
% DuffRHS
%
% Right Hand Side of Duffing equation
%
% Parameters:
%   t      (input)  : current time (not used)
%   u      (input)  : current state u=(x,xp)
%   delta  (input)  : parameter in Duffing equation
%   uprime (output) : function value of right hand side

% -> read input data
x  = u(1);
xp = u(2);

% -> initialize output
uprime = zeros ( size(u) );

% -> evaluate right hand side
uprime(1) = xp;
uprime(2) = x - x*x*x - delta*xp;

