function plotper
% PlotPer
%
% 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 :  Dec 1, 2009
%
% Solution curves of a parameter dependent dynamic system
%
% Parameters
%   none
%
% Example
%   plotper

% -> define initial values, parameter values etc.
mu = 1.538;

tend = 5.0;

x0a = 0.0;
y0a = 2.293331357299218e-001;
x0b = 0.0;
y0b = 2.295615551878201e-001;

ts = 4.0;
ns = 11;
x0vec = linspace ( -0.309, -0.308, ns );
y0vec = -0.56 * ones(size(x0vec));
u0vec = zeros ( 2*ns, 1 );
u0vec(1:2:end) = x0vec;
u0vec(2:2:end) = y0vec;

tolode = 1.0e-12;
options = odeset ( 'AbsTol', tolode, 'RelTol', tolode );

xmin = -0.34;
xmax = -0.28;
ymin = -0.62;
ymax = -0.56;

nquiver = 61;
squiver = 2.0;

% -> time integration
rhs = @(t,u) ( perrhs ( t, u, mu ) );
[ ta, ua ] = ode45 ( rhs, [ 0 tend ], [ x0a y0a ]', options );
[ tb, ub ] = ode45 ( rhs, [ 0 tend ], [ x0b y0b ]', options );

rhs = @(t,uvec) ( mulrhs ( t, uvec, mu ) );
[ tvec, uvec ] = ode45 ( rhs, [ 0 ts ], u0vec, options );

% -> vector field
xx = linspace ( xmin, xmax, nquiver );
yy = linspace ( ymin, ymax, nquiver );
[ XX, YY ] = meshgrid ( xx, yy );

UU = zeros ( size ( XX ) );
VV = zeros ( size ( XX ) );

for ii=1:length(XX(:,1))
    for jj=1:length(XX(1,:))
        uu = perrhs ( 0.0, [ XX(ii,jj) YY(ii,jj) ]', mu );
        UU(ii,jj) = uu(1);
        VV(ii,jj) = uu(2);
    end
end

% -> output
haux = plot ( (xmin+xmax)/2, (ymin+ymax)/2, 'k' );
hold on
quiver ( XX, YY, UU, VV, squiver );
for i=1:ns,
    plot ( uvec(:,2*(i-1)+1), uvec(:,2*(i-1)+2), 'k' )
end
plot ( ua(:,1), ua(:,2), 'g', 'LineWidth', 2 )
plot ( ub(:,1), ub(:,2), 'r', 'LineWidth', 2 )
hold off

set ( haux, 'Visible', 'off' );
axis ( [ xmin xmax ymin ymax ] );
xlabel ( 'x' );
ylabel ( 'y' );
tstr = [ '\mu' sprintf(' = %4.3f',mu) ];
title ( tstr );

end

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

function uprime = perrhs ( t, u, mu )
% PerRHS
%
% Right Hand Side of parameter dependent dynamic system with
% Periodic solution
%
% Parameters:
%   t      (input)  : current time (not used)
%   u      (input)  : current state u=(x,xp)
%   mu     (input)  : system parameter
%   uprime (output) : function value of right hand side

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

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

% -> evaluate right hand side
uprime(1) = mu * ( x - y ) - sin ( x*x + y*y );
uprime(2) = (2+mu)*x + (mu-3)*y - y*y*y;

end

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

function upvec = mulrhs ( t, uvec, mu )
% MulRHS
%
% Right Hand Side of parameter dependent dynamic system with
% periodic solution, version with Multiple right hand sides
%
% Parameters:
%   t      (input)  : current time (not used)
%   uvec   (input)  : current state u=(x,xp), vector valued
%   mu     (input)  : system parameter
%   upvec  (output) : function value of right hand side, vector valued

% -> initialize output
upvec = zeros ( size(uvec) );

% -> evaluate right hand side
for i=1:(length(uvec)/2)
    upvec(2*i+(-1:0)) = perrhs ( t, uvec(2*i+(-1:0)), mu );
end

end
