% graph of a solution to a second-order linear constant
% coefficient equation and its forcing function
%
% as a model of an antihistamine in the bloodstream, we have
%    y'' + (k1 + k2) y' + k2^2 y = (I0/t0) (1 - u(t-t0)),
% or
%    y'' + (k1 + k2) y' + k2^2 y = 
%        (I0/t0) (1 - u(t-t0) + u(t-t1) - u(t-(t1+t0)))
% and so on, for repeated doses.
%
% by Gavin LaRose, 5/10/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 Antihistamine_Graph()
%-------------------------------------------------------------------------------
% configuration variables

% our equation is 
%    y'' + a y' + b y = f(t),
% with initial conditions
%    y(0) = y0, y'(0) = v0
    k1 = 0.128;
    k2 = 0.383;
    k3 = 0.0385;
    
    a = k1 + k2 + k3;
    b = k3*(k1 + k2);
    y0 = 0;
    v0 = 0;
    
% dosing information
    I0 = 25;        % magnitude of dosing
    t0 = 6;         % length of dosing
    t1 = 24;        % time between doses
    num_doses = 2;
    
% t values on which to graph the solutions
    t_range = 0:.02:50;
    
% the forcing function
    function out = f(t)
        out = 0;
        for i = 1:num_doses
            out = out + (I0/t0)*((t>((i-1)*t1)) - (t>((i-1)*t1 + t0)));
        end
        out = k1*out;
    end

% background color
    background_grey = [ 0.6 0.6 0.6 ];

%-------------------------------------------------------------------------------
% figure out the solution
    [t, y] = ode45( @(t,y) [ y(2); -a*y(2) - b*y(1) + f(t) ], ...
                    t_range, [ y0; v0 ] );

    pfig = figure();
    pfaxes = axes( 'Parent', pfig, 'Color', background_grey );
    hold( pfaxes );
    plot( pfaxes, t,y(:,1), '-b', t, f(t), '--k', 'LineWidth', 2 );
    title( [ 'Solution to y'''' + ' num2str(a) ' y'' + ' ...
             num2str(b) ' y = f(t), y(0) = ' num2str(y0) ', y''(0) = ' ...
             num2str(v0) ], 'FontSize', 14 );
    xlabel( 't', 'FontSize', 14 );
    ylabel( 'y', 'FontSize', 14 );
    if num_doses > 1
        legend( 'y', 'f(t)', 'Location', 'NorthWest' );
    else
        legend( 'y', 'f(t)', 'Location', 'SouthEast' );
    end
end

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