% A demonstration to show the solution to the system
%    x1' = -(k01 + k21 + k31) x1 + k12 x2 + k13 x3 + I1
%    x2' = k21 x1 - (k02 + k12) x2
%    x3' = k31 x1 - k13 x3
% with initial condition x1(0) = I01, x2(0) = I02 and x3(0) = I03
%
% The parameters used here are those modeling lead in the body with
% a three-compartment model, so that x1 is the amount of lead in the
% blood, x2 the amount in the tissues and x3 the amount in the bones.
% These are obtained from Borelli, R.L. and C.S. Coleman. Differential
% Equations, a Modeling Perspective, 2nd ed. New York:Wiley, 2004.
%
% by Gavin LaRose, 10/22/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.

function Lead_Euler()
%-------------------------------------------------------------------------------
% configuration variables

% input from the environment
    I1 = 49.3;
    
% initial conditions
    I01 = 0;
    I02 = 0;
    I03 = 0;
    
% transfer coefficients
    k01 = 0.0211;
    k21 = 0.0111;
    k31 = 0.0039;
    k02 = 0.0162;
    k12 = 0.0124;
    k13 = 0.000035;
    
% we show up to the following values of t, then pause; set
%    t_max_second_graph to be the same as t_max_first_graph
%    to omit the second pause
    t_max_first_graph  = 50;
    t_max_second_graph = 1000;
    t_max_final        = 200000;
    
% number of steps for first graph
    num_steps = 50;
% number of steps for second graph
    num_steps_second_graph = 500;
% we increase this for the last graph
    num_steps_final = 5000;
    
% colors
    x1_color = 'r';
    x2_color = 'b';
    x3_color = 'k';
    
    background_grey = [ 0.6 0.6 0.6 ];
    
%-------------------------------------------------------------------------------
% Euler step function

% the transfer matrix and forcing function
    amat = [ -(k01 + k21 + k31) k12          k13; ...
             k21                -(k02 + k12) 0; ...
             k31                0            -k13 ];  
    force = [ I1; 0; 0 ];
    
% a function that generates the Euler approximation for a given number
%    of steps and maximum value of time
    function [ tn, xn ] = euler_approx( n, tmax )
        tn = zeros( 1,n+1 );
        xn = zeros( 3,n+1 );
        xn(1:3,1) = [ I01; I02; I03 ];
        
        h = tmax/n;
        for i=1:n
            tn(i+1) = tn(i) + h;
            xn(1:3,i+1) = xn(1:3,i) + h*(amat*xn(1:3,i) + force);
        end
    end
        
%-------------------------------------------------------------------------------
% initialize graph

    scrsz    = get(0,'ScreenSize');
    scrsz(3) = min(1500, scrsz(3));
    scrsz(4) = min(900, scrsz(4));

    fig1 = figure( 'OuterPosition', ...
                   [50 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2], ...
                   'Name', 'Model Solution' );
    figaxes = axes( 'Parent', fig1, 'Color', background_grey );
    hold on;
    % this ensures we have a toolbar
    set( fig1, 'Toolbar', 'figure' );
    
    xlabel( 't', 'FontSize', 15 );
    ylabel( 'x1,x2,x3', 'FontSize', 15 );
    
    th = title( ['Euler solution 0<=t<=' num2str(t_max_first_graph) ...
                 '; n = ' num2str(num_steps)], 'FontSize', 15 );
    
    
%-------------------------------------------------------------------------------
% Euler solutions

    [tn, xn] = euler_approx( num_steps, t_max_first_graph );
    ph = plot( figaxes, ...
               tn, xn(1,:), [ '.' x1_color ], ...
               tn, xn(2,:), [ '.' x2_color ], ...
               tn, xn(3,:), [ '.' x3_color ], 'MarkerSize', 5 );
    lh = legend( ph, {'x1', 'x2', 'x3'} );
    set( lh, 'FontSize', 15 );
    
    pause on;
    disp( 'Press any key to continue' );
    pause
    
    if t_max_second_graph ~= t_max_first_graph
        % repeat for second max time
        delete( ph );
        delete( lh );
        delete( th );
        
        [tn, xn] = euler_approx( num_steps_second_graph, t_max_second_graph );
        
        th = title( ['Euler solution 0<=t<=' num2str(t_max_second_graph) ...
                     '; n = ' num2str(num_steps_second_graph)], 'FontSize', 15);
    
        ph = plot( tn, xn(1,:), [ '-' x1_color ], ...
                   tn, xn(2,:), [ '-' x2_color ], ...
                   tn, xn(3,:), [ '-' x3_color ], 'LineWidth', 2 );
        lh = legend( ph, {'x1', 'x2', 'x3'}, 'FontSize', 15 );
        
        disp( 'Press any key to continue' );
        pause
    end
    
    % the final graph
    delete( ph );
    delete( lh );
    delete( th );
    
    th = title( ['Euler solution 0<=t<=' num2str(t_max_final) ...
                 '; n = ' num2str(num_steps_final)], 'FontSize', 15 );

    [tn, xn] = euler_approx( num_steps_final, t_max_final );
    
    ph = plot( tn, xn(1,:), [ '-' x1_color ], ...
               tn, xn(2,:), [ '-' x2_color ], ...
               tn, xn(3,:), [ '-' x3_color ], 'LineWidth', 2 );
    lh = legend( ph, {'x1', 'x2', 'x3'}, 'FontSize', 15 );
    
end

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