% euler & improved euler approximation to the system
% x1' = -15 x1 + 4 x2 + 18
% x2' = 3 x1 - 11 x2
%
% by Gavin LaRose, 2/27/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 Lead_Two_Compartment()
%-------------------------------------------------------------------------------
% configuration variables
% the initial condition
    y0 = [ 0; 0 ];
% the step size
    % h = 0.05;
    h  = 0.2;
% the time range on which to find a solution
    t_range = [ 0 1 ];
    
% the transfer matrix and forcing
    coeff_mat = [ -15 4; 3 -11 ];
    forcing   = [ 18; 0 ];
    
% background graph color
    background_grey = [ 0.6 0.6 0.6 ];
    
    if h > 0.1
        constrain_axis1  = 0;
        show_phase_plane = 0;
    else
        constrain_axis1  = 1;
        show_phase_plane = 1;
    end
    % uncomment the following to always omit the phase plane
    % show_phase_plane = 0;
    
    % make this = 0 to omit the pausing between showing the Euler's
    %    and improved Euler methods.
    include_pauses = 1;
    
%-------------------------------------------------------------------------------
% a function giving the RHS

    function out=f(x,y)
        out = coeff_mat*y + forcing;
    end
    
%-------------------------------------------------------------------------------
% calculate the approximations
    nsteps = ceil((t_range(2) - t_range(1))/h);
    
% the Euler approximation
    ye = zeros( 2, nsteps+1 );
    t  = zeros( 1, nsteps+1 );
    
% then ye(1,1:nsteps+1) = the approximations for x1, and ye(2,1:nsteps+1) = x2
    for i=1:nsteps
        t(i+1)      = t(i) + h;
        ye(1:2,i+1) = ye(1:2,i) + h*f(t(i), ye(1:2,i));
    end
    
% the improved Euler approximation
    yi = zeros( 2, nsteps+1 );
    for i=1:nsteps
        k1 = f(t(i), yi(1:2,i));
        u  = yi(1:2,i) + h*k1;
        k2 = f(t(i+1), u);
        yi(1:2,i+1) = yi(1:2,i) + (h/2)*(k1 + k2);
    end
    
%-------------------------------------------------------------------------------
% plot them

    figh  = figure();
    if show_phase_plane
        subp1 = subplot( 1,2,1, 'Parent', figh, 'Color', background_grey );
        subp2 = subplot( 1,2,2, 'Parent', figh, 'Color', background_grey );
    else
        subp1 = axes( 'Parent', figh, 'Color', background_grey );
    end
    title( subp1, 'Approx solutions', ...
           'FontSize', 14 );
    xlabel( subp1, 't', 'FontSize', 14 );
    ylabel( subp1, 'x1, x2', 'FontSize', 14 );
    hold( subp1 );
    
    if show_phase_plane
        title( subp2, 'Phase plane trajectories', ...
               'FontSize', 14 );
        xlabel( subp2, 'x1', 'FontSize', 14 );
        ylabel( subp2, 'x2', 'FontSize', 14 );
        hold( subp2 );
    end
    
    plot( subp1, t,ye(1,1:nsteps+1), '--r', t,ye(2,1:nsteps+1), '--b', ...
          'MarkerSize', 20, 'Marker', '.' );
    if constrain_axis1
        axis( subp1, [ 0 .51 0 1.4 ] );
    else
        % otherwise, make it so that the axes don't resize
        axis(axis);
    end

    if show_phase_plane
        plot( subp2, ye(1,1:nsteps+1),ye(2,1:nsteps+1), '--w', ...
              'MarkerSize', 20, 'Marker', '.' );
    end
    
    if include_pauses
        pause on;
        disp( 'Euler approximation shown.' );
        disp( 'Press any key to plot improved Euler' );
        pause;
    end
    
    plot( subp1, t,yi(1,1:nsteps+1), '-r', 'MarkerFaceColor', 'r', ...
          'MarkerSize', 6, 'Marker', 's' );
    plot( subp1, t,yi(2,1:nsteps+1), '-b', 'MarkerFaceColor', 'b', ...
          'MarkerSize', 6, 'Marker', 's' );

    if show_phase_plane
        plot( subp2, yi(1,1:nsteps+1),yi(2,1:nsteps+1), '-w', ...
              'MarkerFaceColor', 'w', 'MarkerSize', 6, 'Marker', 's' );
    end
end

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