% Modeling the phenomenon of spruce budworm epidemics

% This demo just graphs the solution to the differential equation for
%    initial conditions various
% by Gavin LaRose, 8/9/2012

% (c)2012 Regents of the University of Michigan
% This work is licensed under the
% Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
% To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% or send a letter to
% Creative Commons,
% 444 Castro Street, Suite 900,
% Mountain View, CA, 94041, USA.
%----------------------------------------------------------------------
% configuration variables

% value of q to use
%    these are the values used in the other demos for this section
% interestingQ=[10 9 8 7.3878 7 6];
qval = 10;

% the corresponding colors we use for these solutions
%    these are the values used in the other demos
colorvec={'w','m','y','c','g','b','k'};
col = 'w';

% an initial condition, and a list of initial conditions
ic = 3;
iclist = [ .5, 1, 1.5, 2, 3, 4, 6, 9 ];

% the value of r
r = 0.5;

% where we put arrows on the phase line
phaseLineArrowsLoc = [ .25, 1.3, 3.5, 9 ];

%----------------------------------------------------------------------
% generate the plots (we pause after each one)

scrsz = get(0,'ScreenSize'); % Get the screen size

f1 = figure('Position',[scrsz(3)/4 scrsz(4)/4 scrsz(3)/4 scrsz(4)/2]);
subp1 = subplot( 3, 1, 1, 'Parent', f1, 'Color', [0.5, 0.5, 0.5] );
hold(subp1);
titl = title( ['u'' = ' num2str(r) ' u (1 - u/' num2str(qval) ') - u^2/(1 + u^2)'],...
	      'FontSize', 16 );
subp2 = subplot( 3, 1, 2, 'Parent', f1, 'Color', [0.5, 0.5, 0.5] );
hold(subp2);
subp3 = subplot( 3, 1, 3, 'Parent', f1, 'Color', [0.5, 0.5, 0.5] );
hold(subp3);

% strfmt = ;
pause on

% first put the graph of the solution in the first subplot
for j=1:length(iclist)
    [t,u] = ode45( @(t,u) r.*u.*(1-u/qval)-((u.^2)./(1+u.^2)), [0,200], ...
	iclist(j) );
    if j > length(colorvec)
        colindex = mod(j,length(colorvec));
    else
        colindex = j;
    end
    plot( t, u, [colorvec{colindex} '-'], 'LineWidth', 2, 'Parent', subp1 );
end
xlabel( 't', 'FontSize', 16, 'Parent', subp1 );
ylabel( 'u', 'FontSize', 16, 'Parent', subp1 );
% legend( subp1, ['Solutions with q = ' num2str(qval)], 'FontSize', 16, ...
% 	'Location', 'SouthEast' );
annotation( f1,'textbox', [0.47 0.75 0.14 0.09],...
    'String',{['solutions with various IC']},'FontSize',20,...
    'FitBoxToText','on');

pause

% then plot the corresponding RHS of the ODE in the second subplot
uval = 0:.1:10;
rhsv = r.*uval.*(1 - uval/qval) - ((uval.^2)./(1+uval.^2));
ph=plot( uval, rhsv, 'w-', 'LineWidth', 2, 'Parent', subp2 );
xlabel( 'u', 'FontSize', 16, 'Parent', subp2 );
ylabel( 'u''', 'FontSize', 16, 'Parent', subp2 );
lh = legend( ph, 'RHS of equation', 'Location', 'SouthWest' );
set( lh, 'FontSize', 16 )
plot( uval, zeros(1, length(uval)), [col ':'], 'LineWidth', 2,...
    'Parent', subp2 );

pause

% finally, plot the phase diagram
plot( uval, zeros(1, length(uval)), [col '-'], 'LineWidth', 2,...
    'Parent', subp3 );
for j=1:length(phaseLineArrowsLoc)
    uv = phaseLineArrowsLoc(j);
    if ( r.*uv.*(1 - uv./qval) - ((uv.^2)/(1 + uv.^2)) > 0 )
        plot( uv, 0, [col '>'], 'MarkerSize', 15, 'MarkerFaceColor', ...
	      col, 'Parent', subp3 );
    else
        plot( uv, 0, [col '<'], 'MarkerSize', 15, 'MarkerFaceColor', ...
	      col, 'Parent', subp3 );
    end
end
% add the equilibria
eqn = [ num2str(r), '*x*(1 - x/', num2str(qval), ') - ((x^2)/(1 + x^2))' ];
evals = solve( eqn );
for j=1:length(evals)
    plot( evals(j), 0, [col '.'], 'MarkerSize', 30, 'MarkerFaceColor', ...
	  col, 'Parent', subp3 );
end
annotation( f1,'textbox', [0.4 0.2 0.14 0.09],...
    'String',{['phase diagram']},'FontSize',20,...
    'FitBoxToText','on');

% all done
%----------------------------------------------------------------------
