% Examination of the homogeneous and particular solutions to a second-
% order constant-coefficient linear problem with forcing terms (that
% can be solved with the Method of Undetermined Coefficients).  We
% show two graphs: in the top, the solution to the IVP, y, along with
% the homogeneous and particular solution parts of this; in the bottom,
% the forcing function and the particular solution.
%
% by Gavin LaRose, 10/3/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 MUC_Simple()
%-------------------------------------------------------------------------------
% configuration variables

% we solve a y'' + b y' + c y = f
    a = 1;
    b = 1;
    c = 1;
    
% the function f.  all of the following are good definitions for f
%    f = '2*cos(0.5*t)';
%    f = '2*t*exp(-t/10)';
%    f = '2*exp(-t/10)*cos(2*t)';
    f = '2*cos(0.5*t)';

% initial conditions to use: y(0) = y0, y'(0) = yp0
    y0  = 1;
    yp0 = 0;
    
% the t values on which to graph the solutions
    t_range = [ 0 100 ];
    
% insert pauses between graphing the homogeneous solution, forcing
%    and particular solutions? (0 for no; 1 for yes)
    include_pauses = 1;
    
% color choices
    homog_color  = 'b';
    partic_color = 'g';
    full_color   = 'k';
    force_color  = 'c';
    
    background_grey = [ 0.6 0.6 0.6 ];
    
%-------------------------------------------------------------------------------
% initialization

% we need the vector form of f, too
    fv = regexprep( f, '\*', '.*' );
    fv = regexprep( fv, '\^', '.^' );
    fv = regexprep( fv, '\/', './' );

    scrsz = get(0,'ScreenSize');
    figh = figure('Position',[1 scrsz(4)/4 scrsz(3)/4 scrsz(4)/2]);
    axes( 'Parent', figh, 'Color', background_grey );
    hold on;

% set up subplots
    soln_sub_plot = subplot(2,1,1, 'Parent', figh, 'Color', background_grey );
    forc_sub_plot = subplot(2,1,2, 'Parent', figh, 'Color', background_grey );
    hold( soln_sub_plot );
    hold( forc_sub_plot );
    
    function c=num2charspc(v)
        if v == 1
            c = ' ';
        else
            c = [ num2str(v) ' ' ];
        end
    end
    if b == 0
        term2 = ' ';
    else
        term2 = [ '+ ' num2charspc(b) 'y\prime' ];
    end
    if c == 0
        term3 = ' ';
    else
        term3 = [ '+ ' num2charspc(c) 'y' ];
    end
            
    xlabel( 't', 'FontSize', 16, 'Parent', soln_sub_plot );
    ylabel( 'y', 'FontSize', 16, 'Parent', soln_sub_plot );
    xlabel( 't', 'FontSize', 16, 'Parent', forc_sub_plot );
    ylabel( 'y', 'FontSize', 16, 'Parent', forc_sub_plot );
    title( [ 'Full, homogeneous and particular solutions to' ...
             num2charspc(a) 'y\prime\prime' term2 term3 ' = ' f ], ...
           'FontSize', 16, 'Parent', soln_sub_plot )
    title( [ 'Particular solution and forcing term for' ...
             num2charspc(a) 'y\prime\prime' term2 term3 ' = ' f ], ...
           'FontSize', 16, 'Parent', forc_sub_plot )
    
% the equation we're solving
    eqn = [ num2str(a) '*D2y + ' num2str(b) '*Dy + ' num2str(c) '*y = ' f ];

% define a function for the right-hand side
    function out = rhs(t)
        evalc( [ 'out=' fv ] );
    end
    
%-------------------------------------------------------------------------------
% solve the differential equation and graph the functions

% first we get a function from which we can extricate the particular
%   solution
    % turn off warnings, insofar as we can; we'll test for this later
    warnstate = warning('off','symbolic:mupadmex:MuPADTextWarning');
    ys = dsolve( eqn );
    warnstate = warning('off','symbolic:mupadmex:MuPADTextWarning');
    
% the parameters and variables are in this are
    var_list = symvar( ys );
    
    if length(var_list) > 0
        ygeneral = matlabFunction( ys, 'vars', var_list );

    % for the solution with initial condition, we just use ode45
        [ tvals, ysoln ] = ode45( @(t,y) [ y(2); rhs(t)-(b*y(2)+c*y(1))/a ], ...
                                  t_range, [ y0, yp0 ] );
        
    % all set; do the graphs
        if include_pauses == 1
            pause on
        end
        
    % plot the homogeneous solution, which we define to be the full
    %    solution with the particular solution subtracted off
        plot( tvals, ysoln(:,1) - ygeneral(0,0,tvals), [ '-' homog_color ], ...
              'LineWidth', 2, 'Parent', soln_sub_plot );
        legend1_cells{1} = 'y homogeneous';
        lh1 = legend( soln_sub_plot, legend1_cells );
        
        if include_pauses == 1
            disp( 'Press any key to continue' );
            pause
        end
        
    % next plot the forcing term
        plot( tvals, rhs(tvals), [ '-' force_color ], 'LineWidth', 2, ...
              'Parent', forc_sub_plot );
        legend2_cells{1} = [ 'f(t) = ' f ];
        lh2 = legend( forc_sub_plot, legend2_cells );
        
        if include_pauses == 1
            disp( 'Press any key to continue' );
            pause
        end
        
    % then the full and particular solutions
        plot( tvals, ygeneral(0,0,tvals), [ '-' partic_color ], ...
              tvals, ysoln(:,1), [ '-' full_color ], ...
              'LineWidth', 2, 'Parent', soln_sub_plot );
        legend1_cells{2} = 'y particular';
        legend1_cells{3} = 'full solution y';
        delete(lh1);
        lh1 = legend( soln_sub_plot, legend1_cells );
        
    % and the particular solution in the forcing graph
        plot( tvals, ygeneral(0,0,tvals), [ '-' partic_color ], ...
              'LineWidth', 2, 'Parent', forc_sub_plot );
        legend2_cells{2} = 'y particular';
        delete(lh2);
        lh2 = legend( forc_sub_plot, legend2_cells );
        
    else
        warning(warnstate);
        disp( [ 'Ouch!  Matlab failed to solve the differential equation. ' ...
                 'Is your forcing term one for which we can solve this?' ] );
    end
end

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