function [ sol, T ] = percoll ( u0init, v0init, Tinit, iplot );
% PerColl
%
% This m-file is part of the material of the lecture
%   Scientific Computing and Software Tools
% at the Math Department of Halle University, Germany.
% Summer term 2004
%
% Author :      Prof. Dr. M. Arnold, arnold@mathematik.uni-halle.de
% Version of :  Jun 1, 2004
%
% Computation of a Periodic solution by Matlab's built in collocation method bvp4c
%
% Reference
%   J.D. Murray: Mathematical Biology I, 3rd ed. - Springer - 2002.
%     Section 3.4
%
% Parameter
%   u0init (input)  : initial guess for initial value prey
%   v0init (input)  : initial guess for initial value predator
%   Tinit  (input)  : initial guess for length of one period
%   iplot  (input)  : control flag output
%                       =0 .. no output
%                       =1 .. final output of phase plot
%   sol    (output) : output of boundary value problem solver
%   T      (output) : length of one period
%
% Examples
%
%   [ sol, T ] = percoll ( 0.05, 0.38,  24.0, 1 );
%
%   dbstop in bvp4c at 113;    dbstop in bvp4c at 352;                                                         % set breakpoints
%   [ sol, T ] = percoll ( 0.3, 0.9, 30.0, 1 );                                                                % numerical solution
%   try, myit = it; catch, myit = 0; end;  perplot ( myit, x, y, y(3,end), 50 );   it = myit + 1;    dbcont;   % repeated input in keyboard mode
%   dbclear all                                                                                                % delete breakpoints

% -> error bounds etc.
atol = 1.0e-8;
rtol = 1.0e-4;
nmax = 1.0e4;

% -> get initial guess for bvp4c.m by solving an initial value problem
options = odeset ( 'AbsTol', atol, 'RelTol', rtol );
solode = ode45 ( @wperrhs, [ 0 1 ], [ u0init, v0init, Tinit ]', options, u0init );

Ninit  = 50;
tauini = linspace ( 0, 1, Ninit+1 );
% solinit = bvpinit ( tauini, [ u0init, v0init, Tinit ]' );
solinit.x = tauini;
solinit.y = deval ( solode, tauini );

% -> call the bvp solver
options = bvpset ( 'AbsTol', atol, 'RelTol', rtol, 'NMax', nmax, 'Stat', 'on' );
sol = bvp4c ( @wperrhs, @perbc, solinit, options, u0init );
T = sol.y(3,end);

% -> graphical output (if applicable)
if iplot==1,
    
    tauint = linspace ( 0, 1, 1001 );
    uint   = deval ( sol, tauint, 1 );
    vint   = deval ( sol, tauint, 2 );
    
    subplot ( 2, 1, 1 );
    plot ( uint, vint, 'b', uint(1), vint(1), '+b', uint(end), vint(end), 'ob' );
    axis ( [ 0 0.8 0 0.4 ] );
    xlabel ( 'u' );    ylabel ( 'v' );   
    title ( sprintf ( 'Periodic solution ,  period  T = %8.4f', T ) );

    subplot ( 2, 1, 2 );
    plot ( T*tauint, uint, 'b', T*tauint, vint, 'g' );
    legend ( 'u', 'v' );
    axis ( [ 0 50 0 0.8 ] );
    xlabel ( 't' );
    title ( sprintf ( 'Periodic solution ,  period  T = %8.4f', T ) );
        
end;

% end of PerColl

function Yp = wperrhs ( tau, Y, u0init );
% WPerRHS
%
% Computation of a periodic solution by bvp4c, Wrapper for PerRHS

Yp = perrhs ( tau, Y );

function r = perbc ( Y0, Y1, u0init );
% PerBC
%
% Computation of a periodic solution by bvp4c, Boundary Condition

% -> read input data
u0 = Y0(1);
u1 = Y1(1);
v0 = Y0(2);
v1 = Y1(2);

% -> memory allocation
r = zeros ( size ( Y0 ) );

% -> define output data
r(1) = u0 - u0init;
r(2) = u1 - u0;
r(3) = v1 - v0;

