% 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 + k3) y' + k3(k1 + k2) y = 
%        (I0/t0) (1 - u(t-t0) + u(t-t1) - u(t-(t1+t0)))
% and so on, for repeated doses.
%
% by Gavin LaRose, 5/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 Antihistamine_Doses()
%-------------------------------------------------------------------------------
% configuration variables

% our equation is 
%    y'' + a y' + b y = f(t),
% with initial conditions
%    y(0) = y0, y'(0) = v0
% reasonable values are k1 = 0.128; k2 = 0.383; k3 = 0.0385.  the
%    values here give nice roots to the characteristic equation
%    (roots are r1 = -(k1+k2) and r2 = -k3).
    k1 = 0.12;
    k2 = 0.38;
    k3 = 0.04;
    
    a = k1 + k2 + k3;
    b = k3*(k1 + k2);
    y0 = 0;
    v0 = 0;
    
% dosing information
    I0 = 24;        % magnitude of dosing
    t0 = 6;         % length of dosing
    t1 = 24;        % time between doses
    
% t values on which to graph the solutions
    t_range = 0:.01:250;

% background color
    background_grey = [ 0.6 0.6 0.6 ];

%-------------------------------------------------------------------------------
% figure out the solution
% the forcing function
    function out = f(t)
        out = 0;
        for i = 1:(round(t_range(length(t_range))/t1)+1)
            out = out + (I0/t0)*((t>((i-1)*t1)) - (t>((i-1)*t1 + t0)));
        end
        out = k1*out;
    end

% the base function in our solution
    rvals = roots( [ 1 a b ] );
    r1 = rvals(1);
    r2 = rvals(2);
    function out = g(t)
        out = (k1*I0/t0)*( ( 1/(r1*r2) ) + ...
                           ( 1/(r1*(r1-r2)) )*exp(r1*t) + ...
                           ( 1/(r2*(r2-r1)) )*exp(r2*t));
    end

    t = t_range;
    y = 0*ones(1,length(t));
    for i=1:(round(t(length(t))/t1)+1)
        t_alpha = (i-1)*t1;
        t_beta  = (i-1)*t1 + t0;
        y = y + g(t-t_alpha).*(t>t_alpha) - g(t-t_beta).*(t>t_beta);
    end

    pfig = figure();
    pfaxes = axes( 'Parent', pfig, 'Color', background_grey );
    hold( pfaxes );
    plot( pfaxes, t,y, '-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 );
    legend( 'y', 'f(t)', 'Location', 'NorthWest' );
    axis( [ t_range(1) t_range(length(t_range)) 0 I0 ] );
    axis 'auto y';
end

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