% Examination of the behavior of a second-order linear constant
% coefficient equation, looking at the characteristic polynomial,
% its roots, and the graph of the solution to the differential
% equation.
%
% by Gavin LaRose, 1/2/2013
% parts excerpted from a demonstration by Justin Dunmyre

% (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 Char_Eqn_Simple()
%-------------------------------------------------------------------------------
% configuration variables

% our equation is 
%    m x'' + c x' + k x = 0,
% with initial conditions
%    x(0) = x0, x'(0) = xp0
% good values may be m = 1, c = 0.5, k = 3; and m = 1, c = 3, k = 6
    m = 5;
    c = 22;
    k = 45;
    x0 = 2;
    xp0 = 0;
    
% t values on which to graph the solutions
    t_range = [0 10];
% default y-range on which to graph the solutions
    y_range = [ -5 5 ];
    
% we show numgraphs around the parameters given above, picking values
%    of c that give real distinct, real repeated, and imaginary roots
%    (thus num_graphs has a floor of 3, and should be odd---we add one
%    if it's even)
    num_graphs = 5;
    
% color choices    
    real_colors   = { [0 0 1], [0 .3 1], [0 .5 1], [0 .7 1] };
    repeat_colors = { [1 0 0] };
    imag_colors   = { [0 1 0], [1 1 0], [1 1 .5], [1 1 1] };
    
%-------------------------------------------------------------------------------
% figure out the values to use for c; initially our roots are
    init_roots = roots( [m c k] );
    
    if init_roots(1) == init_roots(2)
        k_start = -c*c/(4*m);
        k_end   = 3*c*c/(4*m);
        k_imag  = k_end;
        real_roots = roots( [m c k_start] );
    else
        k_start = k;
        k_end   = c*c/(2*m) - k;
        if imag(init_roots(1)) == 0
            k_real     = k_start;
            k_imag     = k_end;
            real_roots = init_roots;
            imag_colors = wrev( imag_colors );
        else
            k_real     = k_end;
            k_imag     = k_start;
            real_roots = roots( [m c k_end] );
            real_colors = wrev( real_colors );
        end
    end
    real_roots = sort(real_roots);

% set the number of graphs to produce
    if num_graphs/2 == round(num_graphs/2)
        num_graphs = num_graphs + 1;
    end
    if num_graphs < 3
        num_graphs = 3;
    end
    
% so our increment for k is
    k_incr = (k_end - k_start)/(num_graphs - 1);
    k_values = k_start:k_incr:k_end;
    
%-------------------------------------------------------------------------------
% set up the graph

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

% set up the subplots
    poly_sub_plot = subplot(2,2,1, 'Parent', figh, 'Color',[0.6 0.6 0.6]);
    root_sub_plot = subplot(2,2,2, 'Parent', figh, 'Color',[0.6 0.6 0.6]);
    soln_sub_plot = subplot(2,2,[3 4], 'Parent', figh, 'Color',[0.6 0.6 0.6]);
    hold( poly_sub_plot );
    hold( root_sub_plot );
    hold( soln_sub_plot );
     
% work out the axes for these
    poly_x_range = [ real_roots(1)-1, real_roots(2)+1 ];
    if poly_x_range(1) > 0
        poly_x_range(1) = -1;
    end
    if poly_x_range(2) < 0
        poly_x_range(2) = 1;
    end
    x_mid = (poly_x_range(1) + poly_x_range(2))/2;
    poly_y_min = round(m*x_mid*x_mid + c*x_mid + k_real) - 2;
    tmpval1 = m*poly_x_range(1)*poly_x_range(1) + c*poly_x_range(1) + k_imag;
    tmpval2 = m*poly_x_range(1)*poly_x_range(1) + c*poly_x_range(1) + k_imag;
    poly_y_max   = max([tmpval1 tmpval2]);
    poly_y_range = [ poly_y_min, poly_y_max ];

    rvals1 = sort(roots( [ m c k_start ] ));
    rvals2 = sort(roots( [ m c k_end ] ));
    root_x_range = [ min( [ real(rvals1(1)), real(rvals2(1)) ] ) - 1 ...
                     max( [ real(rvals1(2)), real(rvals2(2)) ] ) + 1 ];
    root_y_range = [ min( [ imag(rvals1(1)), imag(rvals2(1)) ] ) - 1 ...
                     max( [ imag(rvals1(2)), imag(rvals2(2)) ] ) + 1 ];
     
    axis( poly_sub_plot, [ poly_x_range poly_y_range ] );
    axis( root_sub_plot, [ root_x_range root_y_range ] );
    axis( soln_sub_plot, [ t_range y_range ] );
    
% add dashed lines to show axes
    plot( poly_x_range, [0 0], 'w:', 'Parent', poly_sub_plot );
    plot( [0 0], poly_y_range, 'w:', 'Parent', poly_sub_plot );
    plot( root_x_range, [0 0], 'w:', 'Parent', root_sub_plot );
    plot( [0 0], root_y_range, 'w:', 'Parent', root_sub_plot );
    plot( t_range, [0 0], 'w:', 'Parent', soln_sub_plot );
    
% add some titles and axis labels
    poly_title = title( ['Characteristic Eqn: 0 = ', num2str(m), ' r^2 + ' ...
                         num2str(c), ' r + k'], 'FontSize', 16, ...
                        'Parent', poly_sub_plot );
    poly_xlabel = xlabel( 'r', 'FontSize', 16, 'Parent', poly_sub_plot );
    poly_ylabel = ylabel( 'y', 'FontSize', 16, 'Parent', poly_sub_plot );
    root_title = title( 'Roots of Char Eqns', 'FontSize', 16, ...
                        'Parent', root_sub_plot );
    root_xlabel = xlabel( 'Re(r)', 'FontSize', 16, 'Parent', root_sub_plot );
    root_ylabel = ylabel( 'Im(r)', 'FontSize', 16, 'Parent', root_sub_plot );
    soln_title = title( ['Solution to ', num2str(m), 'x\prime\prime + ' ...
                         num2str(c), 'x\prime + k x = 0, x(0) = ', ...
                         num2str(x0), ', x\prime(0) = ', num2str(xp0)], ...
                        'Parent', soln_sub_plot, 'FontSize', 16 );
    soln_xlabel = xlabel( 't', 'FontSize', 16, 'Parent', soln_sub_plot );
    soln_ylabel = ylabel( 'x', 'FontSize', 16, 'Parent', soln_sub_plot );
     
%-------------------------------------------------------------------------------
% populate the figures

    r_vals = poly_x_range(1):.1:poly_x_range(2);
    poly_legend = { };
    root_legend = { };
    soln_legend = { };
    color_index = [0 0 0];
    
    legend_index = 0;
    poly_plot_h = [ ];
    root_plot_h = [ ];
    soln_plot_h = [ ];
    
    % keep track of the range of the solution plot
    y_min = 0;
    y_max = 0;
    for kv=k_values
        legend_index = legend_index + 1;
        
        % figure out the color for the plots
        poly_roots = roots( [m c kv] );
        % ensure that if we're at the case of a repeated root it
        %    doesn't get caught by numerical error
        if legend_index == (num_graphs+1)/2
            poly_roots = [ real(poly_roots(1)), real(poly_roots(1)) ];
        end
        if poly_roots(1) == poly_roots(2)
            i_use = rem(color_index(2), length(repeat_colors))+1;
            use_color = repeat_colors{i_use};
            color_index(2) = color_index(2) + 1;
        elseif imag(poly_roots(1)) == 0
            i_use = rem(color_index(1), length(real_colors))+1;
            use_color = real_colors{i_use};
            color_index(1) = color_index(1) + 1;
        else
            i_use = rem(color_index(3), length(imag_colors))+1;
            use_color = imag_colors{i_use};
            color_index(3) = color_index(3) + 1;
        end
        
        legend_cells{legend_index} = [ 'k = ' num2str(kv) ];
        
        p_h = ...
            plot( r_vals, m*r_vals.*r_vals + c*r_vals + kv, ...
                  'LineStyle', '-', 'Color', use_color, 'LineWidth', 2, ...
                  'Parent', poly_sub_plot );
        poly_plot_h = cat( 2, poly_plot_h, [p_h] );
        r_h = ...
            plot( real(poly_roots), imag(poly_roots), 'Marker', '.', ...
                  'Color', use_color, 'MarkerSize', 20, 'LineStyle', 'none', ...
                  'Parent', root_sub_plot );
        root_plot_h = cat( 2, root_plot_h, [r_h] );
        
        [t_pl, x_pl] = ode45( @(t,x) [ x(2); (-kv*x(1) - c*x(2))/m ], ...
                              t_range, [x0; xp0] );
        s_h = ...
            plot( t_pl, x_pl(:,1), 'LineStyle', '-', 'Color', use_color, ...
              'LineWidth', 2, 'Parent', soln_sub_plot );
        soln_plot_h = cat( 2, soln_plot_h, [s_h] );
        if min(x_pl(:,1)) < y_min
            y_min = min(x_pl(:,1));
        end
        if max(x_pl(:,1)) > y_max
            y_max = max(x_pl(:,1));
        end
        
    end
    if num_graphs > 3
        legend( poly_plot_h, legend_cells, 'Location', 'SouthEastOutside' );
    else
        legend( poly_plot_h, legend_cells, 'Location', 'SouthEast' );
    end
    if num_graphs > 3
        legend( root_plot_h, legend_cells, 'Location', 'SouthEastOutside' );
    else
        legend( root_plot_h, legend_cells, 'Location', 'SouthEast' );
    end
    legend( soln_plot_h, legend_cells, 'Location', 'NorthEast' );
    
    % reset solution axes if needed
    if abs(y_min) < abs(y_range(1))
        y_range(1) = round(y_min - .5);
    end
    if y_max < y_range(2)
        y_range(2) = round(y_max + .5);
    end
    axis( soln_sub_plot, [ t_range, y_range ] );
end

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