function [a_coeff b_coeff c_coeff] = plot_stability_diagrams
%--------------------------------------------------------------------------
% This function plots the Linear Multistep (LM) Implicit-Explicit (ImEx) 
% unconditional stability diagram for a new 1-deltaeter family of ImEx 
% time stepping coefficients.  The new family 
% of coefficients, are deltaeterize by 0 < delta<= 1.  When delta = 1, 
% the coefficients recover the well-known Semi-implicit Backward 
% Differentiation formulas (SBDF). As delta goes to zero, the stability
% region approaches an asymptotic (large) circle.
%
% 09/2016, David Shirokoff
%
% For a given delta, the function also outputs the ImEx time stepping
% coefficients as vectors: (a_coeff, b_coeff, c_coeff). 
% The diagram D, is the set of all point mu in the complex plane, such 
% that the polynomial c(z) - mu b(z) = 0, has every solution |z| < 1.
% The coefficients of the polynomials c(z), b(z) are given by the 
%
% Note: The two variables 'order' and 'delta' are needed to defined
% the ImEx coefficiets.  
%--------------------------------------------------------------------------

    %The new coefficients are defined for orders 1-5, and delta parameter 
    %0 < delta <= 1. Modify these variables to plot diagrams for schemes 
    %with different orders.
    order = 2;  
    delta = 1;  
    
    plot_asymptotic_circle = 0;     %1-plots the asymptotic circle
                                    %0-turns circle off.
                                    
    %The angle of the complex number z0 is used in the boundary local 
    %method.  z0 describes the starting and ending
    %points of the parameterized boundary D.
    z0 = ((2-delta) - 2*(1-delta)*cos(pi/order)*exp(i*pi/order))./(2-delta - 2*cos(pi/order)*exp(i*pi/order));
    ang_curve = angle(z0);
    [a_coeff b_coeff c_coeff a b c] = GetImExCoefficients(order, delta, ang_curve);
    
    LW = 'linewidth'; lw = 3; MS = 'markersize'; ms = 8;
    hold on;
    v = c./b;
    plot(real(v), imag(v), 'Color', 'k', MS, ms, LW, 5);
    grid on, box on, axis square
    set(gca,'FontSize',30, LW, 3)
    
    title( strcat(strcat('Stability diagram, Order= ', num2str(order)), ', \delta= ', num2str(delta)))
    fill(real(v), imag(v), 'b')
    label2 = 'D';
    label1 = '\partial D';
    if (plot_asymptotic_circle == 1)
        th = 0:0.1:2*pi;
        z = exp(i*th);
        w = 1/(order*delta)*(z - 1) + (order + 1)/(2*order);
        plot(real(w), imag(w), 'o', 'Color', 'k', LW, 2)
        label3 = 'Asy. Circle';
        legend(label1,label2, label3)
    else
        label1 = '\partial D';
        legend(label1, label2)
    end    
end

%------------------------------------------------------------------
% Function takes the value of order and delta
% and outputs the ImEx coefficients and polynomials evaluated around
% the unit circle z = e^(i\theta).
%
function [a_coeff b_coeff c_coeff a b c] = GetImExCoefficients(order, delta, ang_curve)
    %Note: a, b, c variables are vectors of the polynomials a(z), b(z) and
    %c(z) evaluated at values of z on a subset of the unit circle.
    th = ang_curve:0.001:(2*pi - ang_curve);
    z = exp(i*th);
    switch order        
        case 1
            a = delta*z - delta; 
            c = delta + z - 1;  
            b = delta;
            
            a_coeff = [delta -delta];
            b_coeff = [0 delta];
            c_coeff = [1 delta-1];

        case 2
            a = (2*delta - delta^2/2)*z.^2 + (2*delta^2 - 4*delta)*z + 2*delta - (3*delta^2)/2;
            c = z*(2*delta - 2) + (delta - 1)^2 + z.^2;
            b = (delta - 1)^2 + 2*delta*z - 1;
            
            a_coeff = [(2*delta - delta^2/2) (2*delta^2 - 4*delta) (2*delta-(3*delta^2)/2)];
            b_coeff = [0 2*delta (delta - 1)^2 - 1];
            c_coeff = [1 (2*delta - 2) (delta - 1)^2];
            
        case 3
            a = z*(3*delta^3 - (21*delta^2)/2 + 9*delta) - 3*delta + z.^3*(delta^3/3 - (3*delta^2)/2 + 3*delta) - z.^2*((3*delta^3)/2 - (15*delta^2)/2 + 9*delta) + (9*delta^2)/2 - (11*delta^3)/6;
            c = (delta - 1)^3 + z.^2*(3*delta - 3) + z*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + z.^3;
            b = (delta - 1)^3 + 3*delta*z.^2 + z*((delta - 1)^2 + (2*delta - 2)*(delta - 1) - 3) + 1;
            
            a_coeff =[(delta^3/3 - (3*delta^2)/2 + 3*delta), -((3*delta^3)/2 - (15*delta^2)/2 + 9*delta), (3*delta^3 - (21*delta^2)/2 + 9*delta), (-3*delta + (9*delta^2)/2 - (11*delta^3)/6) ]
            b_coeff =[ 0, 3*delta,  ((delta - 1)^2 + (2*delta - 2)*(delta - 1) - 3), (delta - 1)^3 + 1]
            c_coeff = [1, (3*delta - 3), ((delta - 1)^2 + (2*delta - 2)*(delta - 1)), (delta - 1)^3]
            
        case 4 
            a = 4*delta - z*(- 4*delta^4 + (58*delta^3)/3 - 30*delta^2 + 16*delta) - 9*delta^2 + (22*delta^3)/3 - (25*delta^4)/12 + z.^4*(- delta^4/4 + (4*delta^3)/3 - 3*delta^2 + 4*delta) - z.^3*(- (4*delta^4)/3 + (22*delta^3)/3 - 18*delta^2 + 16*delta) + z.^2*(- 3*delta^4 + 18*delta^3 - 36*delta^2 + 24*delta);
            c = (delta - 1)^4 + z.^3*(4*delta - 4) + z.^2*((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1)) + z*((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3) + z.^4; 
            b = (delta - 1)^4 + z.^2*((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1) - 6) + z*((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3 + 4) + 4*delta*z.^3 - 1;
            
            a_coeff = [(- delta^4/4 + (4*delta^3)/3 - 3*delta^2 + 4*delta), -(-(4*delta^4)/3 + (22*delta^3)/3 - 18*delta^2 + 16*delta), (- 3*delta^4 + 18*delta^3 - 36*delta^2 + 24*delta), -(- 4*delta^4 + (58*delta^3)/3 - 30*delta^2 + 16*delta), 4*delta - 9*delta^2 + (22*delta^3)/3 - (25*delta^4)/12 ]
            b_coeff = [ 0, 4*delta, ((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1) - 6),  ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3 + 4), (delta - 1)^4  - 1];
            c_coeff = [1, (4*delta - 4), ((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1)), ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3), (delta - 1)^4]; 
            
        case 5
            a = z.^5*(delta^5/5 - (5*delta^4)/4 + (10*delta^3)/3 - 5*delta^2 + 5*delta) - 5*delta - z.^4*((5*delta^5)/4 - (95*delta^4)/12 + (65*delta^3)/3 - 35*delta^2 + 25*delta) + z.^3*((10*delta^5)/3 - (65*delta^4)/3 + (190*delta^3)/3 - 90*delta^2 + 50*delta) - z.^2*(5*delta^5 - 35*delta^4 + (280*delta^3)/3 - 110*delta^2 + 50*delta) + z*(5*delta^5 - (365*delta^4)/12 + (200*delta^3)/3 - 65*delta^2 + 25*delta) + 15*delta^2 - (55*delta^3)/3 + (125*delta^4)/12 - (137*delta^5)/60;
            c = (delta - 1)^5 + z*((delta - 1)^4 + ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3)*(delta - 1)) + z.^2*((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3 + (delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1))) + z.^4*(5*delta - 5) + z.^5 + z.^3*((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1) + (4*delta - 4)*(delta - 1)); 
            b =(delta - 1)^5 + z*((delta - 1)^4 + ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3)*(delta - 1) - 5) + z.^2*((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3 + (delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1)) + 10) + 5*delta*z.^4 + z.^3*((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1) + (4*delta - 4)*(delta - 1) - 10) + 1;
 
            a_coeff =[(delta^5/5 - (5*delta^4)/4 + (10*delta^3)/3 - 5*delta^2 + 5*delta), -((5*delta^5)/4 - (95*delta^4)/12 + (65*delta^3)/3 - 35*delta^2 + 25*delta), ((10*delta^5)/3 - (65*delta^4)/3 + (190*delta^3)/3 - 90*delta^2 + 50*delta), -(5*delta^5 - 35*delta^4 + (280*delta^3)/3 - 110*delta^2 + 50*delta), (5*delta^5 - (365*delta^4)/12 + (200*delta^3)/3 - 65*delta^2 + 25*delta), 15*delta^2 - (55*delta^3)/3 + (125*delta^4)/12 - (137*delta^5)/60 - 5*delta];
            b_coeff = [0, 5*delta, ((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1) + (4*delta - 4)*(delta - 1) - 10),  ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3 + (delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1)) + 10), ((delta - 1)^4 + ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3)*(delta - 1) - 5), 1 + (delta - 1)^5];
            c_coeff = [1, (5*delta - 5), ((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1) + (4*delta - 4)*(delta - 1)), ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3 + (delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1))), ((delta - 1)^4 + ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3)*(delta - 1)), (delta - 1)^5];

    end
end