% Examination of the behavior of a second- or third-order linear
% constant coefficient equation, looking at the characteristic
% equation, its roots, and the graph of the solution to the
% differential equation.  This demo animates the variation of one
% parameter to show the effect of this on the characteristic equation,
% the roots, and the solution.
%
% This solves the differential equation
%    a x''' + b x'' + c x' + d x = 0,
% with initial conditions [ x(0); x'(0); x''(0) ] = initCond
%
% by Justin Dunmyre
% lightly edited by Gavin LaRose 9/19/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 Char_Eqn_GUI()
%-------------------------------------------------------------------------------
% configuration variables

% these give the initial valus of the parameters a, b, c and d as shown
%    in the GUI control box.
    defaulta = 0;
    defaultb = 1;
    defaultc = 0;
    defaultd = 1;

% the initial conditions, x(0), x'(0) and x''(0)
    initCond = [1; 1; 0];

% the parameter to vary in the animation
    parToVary = 'd';

% the minimum and maximum values of the varying parameters, and the 
%    number of steps to take in the animation.
    parmin = -3;
    parmax = 3;
    numpar = 80;
    
  
% the end of the range of t-values considered in the solution of the
%    differential equation
    tend = 30;
    
% the color to use when graphing the axes (where x = y = 0)
    zeroColor = 'c';
% the color to use for the graphs
    solColor  = 'k';
    
% the background color
    bggray=[.6 .6 .6];
% the range of values on which to graph the characteristic equation
    rRange=linspace(-10,10,10000);

% the range and domain used when graphing the characteristic equation
    charAxesX=[rRange(1) rRange(end)];
    charAxesY=[-5 5];
    
% the range and domain used when graphing the solution to the
%    differential equation
    solAxesX=[0 tend];
    solAxesY=[-5 5];

%-------------------------------------------------------------------------------
% initialization

    a = defaulta;
    b = defaultb;
    c = defaultc;
    d = defaultd;
    
    parRange = linspace(parmin,parmax,numpar);
    
    togglebutton = 1;

    polySubplot  = 0;
    solSubplot   = 0;
    rootsSubplot = 0;
    figh         = 0;
    controls     = 0;
    initRootPlot = 0;
    initSolPlot  = 0;
    initCharPlot = 0;
    myplot       = 0;
    mysolPlot    = 0;
    rootplot     = 0;
    scrsz        = 0;
    
    ta = 0;
    tb = 0;
    tc = 0;
    td = 0;
    
    drawWindow();
    
    hold on;
    
    pauseFlag = 0;
    
    setupUI();
    
%-------------------------------------------------------------------------------
% the main processing
    
function mainProcess()    

    figure(figh)
    title(solSubplot,formODEstring(a,b,c,d),'FontSize',15);
    set(initRootPlot,'Visible','off');
    set(initCharPlot,'Visible','off');
    set(initSolPlot,'Visible','off');
    set(myplot(end),'Visible','off');
    set(mysolPlot(end),'Visible','off');
    set(rootplot(end),'Visible','off');
    
    
    getframe();
    y=zeros(length(parRange),length(rRange));
        
    for j=1:length(parRange)
        parVal=parRange(j);
        evalc([parToVary '=' num2str(parVal)]); %parToVary=parVal, used this way to let parToVary be any letter
        
        y(j,:)=charPolynomial(rRange);
        
        [tempt,tempsol]=ode45(@myODE,[0 tend],initCond);
        
        t{j}=tempt;
        sol{j}=tempsol(:,1);
        if rem(j,10) == 0
            disp(['Solving ODE: ' num2str(j) ' of ' ...
                  num2str(length(parRange))]);
        end
    end
    
    parVal=parRange(1);
    myplot(1)=plot(rRange,y(j,:),[solColor '-'],'LineWidth',2,'Parent',polySubplot);

    mysolPlot(1)=plot(t{1},sol{1},[solColor '-'],'LineWidth',2,'Parent',solSubplot);
    
    title(solSubplot,formODEstring(a,b,c,d),'FontSize',15);
    
    myRoots=roots([a b c d]);

    rootplot(1)=plot(real(myRoots),imag(myRoots),[solColor '.'],'MarkerSize',20,'Parent',rootsSubplot);
    
    getframe();
        
    for j=2:length(parRange)

        set(myplot(j-1),'Visible','off');
        set(mysolPlot(j-1),'Visible','off');
        set(rootplot(j-1),'Visible','off');
        parVal=parRange(j);
        evalc([parToVary '=' num2str(parVal)]);
        
        myplot(j)=plot(rRange,y(j,:),[solColor '-'],'LineWidth',2,'Parent',polySubplot);
        legend(polySubplot,[parToVary ' = ' num2str(parVal)]);

        
        mysolPlot(j)=plot(t{j},sol{j},[solColor '-'],'LineWidth',2,'Parent',solSubplot);
        title(solSubplot,formODEstring(a,b,c,d),'FontSize',15);    


        myRoots=roots([a b c d]);
        rootplot(j)=plot(real(myRoots),imag(myRoots),[solColor '.'],'MarkerSize',20,'Parent',rootsSubplot);

        getframe();
        
    end
end

%-------------------------------------------------------------------------------
% functions for the characteristic polynomial and differential equation       
       
function out=charPolynomial(x)
    out=a*x.^3+b*x.^2+c*x+d;
end

function LHS=myODE(t,xvec)
%ax'''+bx''+cx'+dx=0

    x1=xvec(1);
    x2=xvec(2);
    x3=xvec(3);
       
    if(a~=0)
        x1prime=x2;
        x2prime=x3;
        x3prime=(-b*x3-c*x2-d*x1)/a;
    elseif(a==0)
        x1prime=x2;
        x2prime=(-c*x2-d*x1)/b;
        x3prime=0;
    end
       
    LHS=[x1prime;x2prime;x3prime];
end
   
%-------------------------------------------------------------------------------
% Return a string for the differential equation

function stringy=formODEstring(odea,odeb,odec,oded)
    stringy=[];
    if(odea~=0)
        stringy=[num2str(odea) ' x\prime\prime\prime'];
    end
    if(odeb~=0)
        if(odea~=0)
            stringy=[stringy ' + ' num2str(odeb) 'x\prime\prime'];
        else
            stringy=[num2str(odeb) 'x\prime\prime'];
        end
    end
    if(c~=0)
        if(odea~=0||odeb~=0)
            stringy=[stringy ' + ' num2str(odec) 'x\prime'];
        else
            stringy=[num2str(odec) 'x\prime'];
        end
    end
    if(oded~=0)
        stringy=[stringy ' + ' num2str(oded) 'x'];
    end
    stringy=[stringy '=0'];
end

%-------------------------------------------------------------------------------
% Control of radial buttons a b c d
%    when a radial button is selected, it hides the corresponding
%    parameter edit area

function uibuttongroup1_SelectionChangeFcn(hObject,eventdata)
    switch get(eventdata.NewValue,'Tag') % Get Tag of selected object.
      case 'a'
        set(ta,'Style','text','ForegroundColor',get(ta,'BackgroundColor'));
        set(tb,'Style','edit','ForegroundColor','k');
        set(tc,'Style','edit','ForegroundColor','k');
        set(td,'Style','edit','ForegroundColor','k');
        parToVary='a';
        % Code for when radiobutton1 is selected.
      case 'b'
        set(ta,'Style','edit','ForegroundColor','k');
        set(tb,'Style','text','ForegroundColor',get(ta,'BackgroundColor'));
        set(tc,'Style','edit','ForegroundColor','k');
        set(td,'Style','edit','ForegroundColor','k');
        parToVary='b';
        % Code for when radiobutton2 is selected.
      case 'c'
        set(ta,'Style','edit','ForegroundColor','k');
        set(tb,'Style','edit','ForegroundColor','k');
        set(tc,'Style','text','ForegroundColor',get(ta,'BackgroundColor'));
        set(td,'Style','edit','ForegroundColor','k');
        parToVary='c';
        % Code for when togglebutton1 is selected.
      case 'd'
        set(ta,'Style','edit','ForegroundColor','k');
        set(tb,'Style','edit','ForegroundColor','k');
        set(tc,'Style','edit','ForegroundColor','k');
        set(td,'Style','text','ForegroundColor',get(ta,'BackgroundColor'));
        parToVary='d';
        % Code for when togglebutton2 is selected.
        % Continue with more cases as necessary.
      otherwise
    end
end

%-------------------------------------------------------------------------------
% start button

function startButtonCall(hObject,eventdata)
            
    set(initRootPlot,'Visible','off');
    set(initCharPlot,'Visible','off');
    set(initSolPlot,'Visible','off');
    set(myplot(end),'Visible','off');
    set(mysolPlot(end),'Visible','off');
    set(rootplot(end),'Visible','off');  
    
    mainProcess();
    
end

%-------------------------------------------------------------------------------
% pause / continue button

function togglebuttonCall(hObject,eventdata)
    if(pauseFlag==0)
        set(togglebutton,'String','Continue');
        pauseFlag=1;
        uiwait;
    else
        pauseFlag=0;
        set(togglebutton,'String','Pause');
        uiresume;
        figure(figh)
    end
end

%-------------------------------------------------------------------------------
% reset button

function resetButtonCall(hObject,eventdata)

    uiresume;
     
    a = defaulta;
    b = defaultb;
    c = defaultc;
    d = defaultd;
    
    set(initRootPlot,'Visible','off');
    set(initCharPlot,'Visible','off');
    set(initSolPlot,'Visible','off');
    set(myplot(end),'Visible','off');
    set(mysolPlot(end),'Visible','off');
    set(rootplot(end),'Visible','off');  
    
    pauseFlag=0;
    parToVary='d';
    setupUI();
end

%-------------------------------------------------------------------------------
% update appropriate variables when their text box changes

function t1cb(hObject,eventdata)
    a=str2double(get(hObject,'String'));
end
function t2cb(hObject,eventdata)
    b=str2double(get(hObject,'String'));
end
function t3cb(hObject,eventdata)
    c=str2double(get(hObject,'String'));
end
function t4cb(hObject,eventdata)
    d=str2double(get(hObject,'String'));
end

function parminCB(hObject,eventdata)
    parmin=str2double(get(hObject,'String')); 
    parRange=linspace(parmin,parmax,numpar);
end
function parmaxCB(hObject,eventdata)
    parmax=str2double(get(hObject,'String')); 
    parRange=linspace(parmin,parmax,numpar);
end

%-------------------------------------------------------------------------------
% initialize graphs button

function initGraphCallback(hObject,eventdata)
    %hide previous plots
        set(initRootPlot,'Visible','off');
        set(initCharPlot,'Visible','off');
        set(initSolPlot,'Visible','off');
        set(myplot(end),'Visible','off');
        set(mysolPlot(end),'Visible','off');
        set(rootplot(end),'Visible','off');   
    
    %generate information for new plots
        evalc([parToVary '=' num2str(parmin)]);
        cPoly=charPolynomial(rRange);
        myRoots=roots([a b c d]);
        [myt,mysol]=ode45(@myODE,[0 tend],initCond);
    %make new plots
        initCharPlot=plot(rRange,cPoly,[solColor '-'],'LineWidth',2,'Parent',polySubplot);
        initRootPlot=plot(real(myRoots),imag(myRoots),[solColor '.'],'MarkerSize',20,'Parent',rootsSubplot);        
        legend(polySubplot,[parToVary ' = ' num2str(parmin)]);
        initSolPlot=plot(myt,mysol(:,1),[solColor '-'],'LineWidth',2,'Parent',solSubplot);
        title(solSubplot,formODEstring(a,b,c,d),'FontSize',15);    
end

%-------------------------------------------------------------------------------
% set up the UI

function setupUI()
    
    bgroup = uibuttongroup('visible','off','Position',[0 0 1 1],'Parent',controls);
    % Create three radio buttons in the button group.
    togglebutton = uicontrol('Parent',bgroup,'Style','PushButton','Min',0,'Max',1,'Parent',bgroup,...
                             'Position',[50 0 100 30],'String','Pause/Continue','Callback',@togglebuttonCall);
        
    endbutton= uicontrol('Parent',bgroup,'Style','PushButton','Parent',bgroup,...
                         'Position',[50 480 100 30],'String','Reset','Callback',@resetButtonCall);
        
        
    pushbutton = uicontrol('Parent',bgroup,'Style','PushButton','Parent',bgroup,...
                           'Position',[5 515 90 30],'String','Start','Callback',@startButtonCall);
    initButton = uicontrol('Parent',bgroup,'Style','PushButton','pos',[105 515 90 30],'String','Initialize Graphs','Callback',@initGraphCallback);
    
        
   
    u0 = uicontrol('Style','Radio','String','a',...
                   'pos',[10 370 100 30],'parent',bgroup,'HandleVisibility','off','Tag','a');
    u1 = uicontrol('Style','Radio','String','b',...
                   'pos',[50 370 100 30],'parent',bgroup,'HandleVisibility','off','Tag','b');
    u2 = uicontrol('Style','Radio','String','c',...
                   'pos',[90 370 100 30],'parent',bgroup,'HandleVisibility','off','Tag','c');
    u3 = uicontrol('Style','Radio','String','d',...
                   'pos',[130 370 100 30],'parent',bgroup,'HandleVisibility','off','Tag','d');
        
        
    uicontrol('Style','text','parent',bgroup,'pos',[10 250 30 30],'String','a','FontSize',15);
    uicontrol('Style','text','parent',bgroup,'pos',[10 200 30 30],'String','b','FontSize',15);
    uicontrol('Style','text','parent',bgroup,'pos',[10 150 30 30],'String','c','FontSize',15);
    uicontrol('Style','text','parent',bgroup,'pos',[10 100 30 30],'String','d','FontSize',15);
    ta=uicontrol('Style','edit','parent',bgroup,'pos',[50 255 30 30],'Callback',@t1cb,'String',num2str(a));
    tb=uicontrol('Style','edit','parent',bgroup,'pos',[50 205 30 30],'Callback',@t2cb,'String',num2str(b));
    tc=uicontrol('Style','edit','parent',bgroup,'pos',[50 155 30 30],'Callback',@t3cb,'String',num2str(c));
    td=uicontrol('Style','edit','parent',bgroup,'pos',[50 105 30 30],'Callback',@t4cb,'String',num2str(d));
        
    c1= uicontrol('Style','text','parent',bgroup,'HandleVisibility','off',...
                  'pos',[0 430 180 30],'String','ax''''''+bx''''+cx''+dx=0','FontSize',15);
    c2= uicontrol('Style','text','parent',bgroup,'HandleVisibility','off',...
                  'pos',[0 400 180 30],'String','Select a,b,c,d to vary that parameter (default d)','FontSize',10);
    uicontrol('Style','text','parent',bgroup,'pos',[40 330 70 30],'String','between','FontSize',15);
    
    uicontrol('Style','edit','parent',bgroup,'pos',[20 305 30 30],'Callback',@parminCB,'String',num2str(parmin));
    uicontrol('Style','text','parent',bgroup,'pos',[50 300 50 30],'String','and','FontSize',15);
    uicontrol('Style','edit','parent',bgroup,'pos',[100 305 30 30],'Callback',@parmaxCB,'String',num2str(parmax));
    uicontrol('Style','text','parent',bgroup,'pos',[5 60 190 40],'String',['x(0)=' num2str(initCond(1)) ' x''(0)=' num2str(initCond(2)) ' and if a =/= 0, then x''''(0)=' num2str(initCond(3))],'FontSize',15);
    % Initialize some button group properties. 
        
    switch(parToVary)
      case 'a'
        set(ta,'Style','text','ForegroundColor',get(ta,'BackgroundColor'));
        selected=u0;
      case 'b'
        set(tb,'Style','text','ForegroundColor',get(ta,'BackgroundColor'));
        selected=u1;
      case 'c'
        set(tc,'Style','text','ForegroundColor',get(ta,'BackgroundColor'));                
        selected=u2;
      case 'd'
        set(td,'Style','text','ForegroundColor',get(ta,'BackgroundColor'));
        selected=u3;
      otherwise
        selected=[];
    end
    set(bgroup,'SelectionChangeFcn',@uibuttongroup1_SelectionChangeFcn);
    set(bgroup,'SelectedObject',[selected]);  % No selection
    set(bgroup,'Visible','on');   
           
end

%-------------------------------------------------------------------------------
% create figures

function drawWindow()
    scrsz = get(0,'ScreenSize'); %Get the screen size
    pause on %just to be sure, enable pausing

    controls=figure('OuterPosition',[1 scrsz(4)/5 205 scrsz(4)/1.5]);
    set(controls,'MenuBar','none','Name','Controls');
    figh=figure('OuterPosition',[200 scrsz(4)/3 scrsz(3)/3 scrsz(4)/2]); %Create a new figure window
    axes('Parent',figh,'Color',bggray); %Change the color


    polySubplot=subplot(2,2,1,'Parent',figh,'Color',bggray);
    solSubplot=subplot(2,2,[3 4],'Parent',figh,'Color',bggray);
    rootsSubplot=subplot(2,2,2,'Parent',figh,'Color',bggray);

    hold(polySubplot);
    hold(solSubplot);
    hold(rootsSubplot);

    axis(polySubplot,[charAxesX charAxesY]);
    axis(solSubplot,[solAxesX solAxesY]);
    axis(rootsSubplot,[-2 2 -2 2]); 
    title(polySubplot,'Characteristic Polynomial','FontSize',15);
    xlabel(polySubplot,'r','FontSize',15);

    p=plot(linspace(-2,2,100),zeros(1,100),[zeroColor ':'],'LineWidth',2,'Parent',rootsSubplot);
    set(get(get(p,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','off'); 
    p=plot(zeros(1,100),linspace(-2,2,100),[zeroColor ':'],'LineWidth',2,'Parent',rootsSubplot);
        set(get(get(p,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','off');
    
    p=plot(rRange,zeros(1,length(rRange)),[zeroColor ':'],'LineWidth',2,'Parent',polySubplot,'DisplayName','');
        set(get(get(p,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','off');
    p=plot(zeros(1,length(rRange)),linspace(-5,5,length(rRange)),[zeroColor ':'],'LineWidth',2,'Parent',polySubplot,'DisplayName','');
        set(get(get(p,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','off');
    xlabel(rootsSubplot,'Re','FontSize',15);
    ylabel(rootsSubplot,'Im','FontSize',15);
    title(rootsSubplot,'Roots of Char. Poly.','FontSize',15);
end

end

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