% Plot solutions to the logistic and spruce budworm model for
% different initial conditions, showing how the behaviors of 
% the two are similar and different.  We plot
%    P' = k P (M - P)
% and 
%    P' = k P (M - P) - P^2/(1 + P^2)
% 
% by Gavin LaRose, 1/15/2013
%
% (c)2013 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.

function Logistic_Budworm()
%----------------------------------------------------------------------
% configuration variables
    k = 0.05;
    M = 10;
    
    t_range = [ 0 50 ];
% the first initial condition is shown, then the others added after
%    a pause
    p0_vals = [ 3 1 5 8 10 13 15 ];
    
    background_grey = [ 0.6 0.6 0.6 ];
%----------------------------------------------------------------------
% generate some numerical solutions
% first show the solution for the first initial condition
    tvals = t_range(1):.05:t_range(2);
    [t, ls1] = ode45( @(t,p) k*p.*(M-p), tvals, p0_vals(1) );
    
    figh = figure();
    axes( 'Parent', figh, 'Color', background_grey);
    axis( [ 0 50 0 15 ] );
    xlab = xlabel( 't', 'FontSize', 16 );
    ylab = ylabel( 'P', 'FontSize', 16 );
    titl = title( ['Solutions to P'' = k P (M - P)'],...
            'FontSize', 18 );
    hold on;
    pause on;
    
    plot( t, ls1, '-b', 'LineWidth', 2 );
    disp( 'Press any key to show other solutions' );
    pause;
    
    for i = 2:length(p0_vals)
        [t, ls] = ode45( @(t,p) k*p.*(M - p), tvals, p0_vals(i) );
        plot( t, ls, '--k', 'LineWidth', 2 );
    end
    
    disp( 'Press any key to show solutions to budworm equation' );
    pause;
    
    [t, bs1] = ode45( @(t,p) k*p.*(M-p) - p.^2./(1+p.^2), tvals, ...
        p0_vals(1) );
    
    figg = figure();
    ah = axes( 'Parent', figg, 'Color', background_grey);
    axis( [ 0 50 0 15 ] );
    xlab = xlabel( 't', 'FontSize', 16 );
    ylab = ylabel( 'P', 'FontSize', 16 );
    titl = title( ['Solutions to P'' = k P (M - P) - P^2/(1 + P^2)'],...
            'FontSize', 18 );
    hold on;
    pause on;
    
    plot( ah, t, bs1, '-b', 'LineWidth', 2 );
    disp( 'Press any key to show other solutions' );
    pause;
    
    for i = 2:length(p0_vals)
        [t, bs] = ode45( @(t,p) k*p.*(M - p) - p.^2./(1 + p.^2), ...
            tvals, p0_vals(i) );
        plot( ah, t, bs, '--k', 'LineWidth', 2 );
    end
end
