function temple5044_stability_region_imex_rk
%TEMPLE5044_STABILITY_REGION_IMEX_RK
%   Plots the stability regions for a number of
%   ImEx Runge-Kutta schemes (with DIRK), for the
%   test problem u'=alpha*u+i*beta*u.
%
%   (C) 2019/04/02 by Benjamin Seibold
%       http://www.math.temple.edu/~seibold/

%========================================================================
% Define RK schemes
%========================================================================
name{1} = 'forward backward Euler (1,1,1)';
A{1} = [1];
b{1} = [1];
Ahat{1} = [0 0;1 0];
bhat{1} = [1;0];

name{2} = 'forward backward Euler (1,2,1)';
A{2} = [1];
b{2} = [1];
Ahat{2} = [0 0;1 0];
bhat{2} = [0;1];

name{3} = 'ImEx midpoint (1,2,2)';
A{3} = [1/2];
b{3} = [1];
Ahat{3} = [0 0;1/2 0];
bhat{3} = [0;1];

name{4} = 'Third order (2,3,3)';
gamma = (3+sqrt(3))/6;
A{4} = [gamma 0;1-2*gamma gamma];
b{4} = [1/2;1/2];
Ahat{4} = [0 0 0;gamma 0 0;gamma-1 2*(1-gamma) 0];
bhat{4} = [0;1/2;1/2];

name{5} = 'L-stable (2,2,2)';
gamma = (2-sqrt(2))/2; delta = 1-1/(2*gamma);
A{5} = [gamma 0;1-gamma gamma];
b{5} = A{5}(end,:)';
Ahat{5} = [0 0 0;gamma 0 0;delta 1-delta 0];
bhat{5} = Ahat{5}(end,:)';

name{6} = 'L-stable (4,4,3)';
A{6} = [1/2 0 0 0;1/6 1/2 0 0;-1/2 1/2 1/2 0;3/2 -3/2 1/2 1/2];
b{6} = A{6}(end,:)';
Ahat{6} = [0 0 0 0 0;1/2 0 0 0 0;...
    11/18 1/18 0 0 0;5/6 -5/6 1/2 0 0;1/4 7/4 3/4 -7/4 0];
bhat{6} = Ahat{6}(end,:)';

%========================================================================
% Plot stability regions
%========================================================================
ax = 9; % size of stability region plot
n = length(name); % number of schemes
nc = ceil(sqrt(n)); % number of subplot columns
nr = ceil(n/nc); % number of subplot rows
clf
for j = 1:length(name) % loop over schemes
    subplot(nr,nc,j)
    plot_stability_region(A{j},b{j},Ahat{j},bhat{j},ax) % plot region
    hold on
    plot([-ax,ax],[0,0],'k:',[0,0],[-ax,ax],'k:')
    hold off
    title(name{j})
end

%========================================================================
% Functions
%========================================================================
function plot_stability_region(A,b,Ahat,bhat,ax)
if nargin<5, ax = 1; end % size of stability region plot
s = length(b); % number of implicit stages
zr = linspace(-1,1,200)*ax; zi = zr;
[ZR,ZI] = meshgrid(zr,zi);
R = ZR*0;
for jr = 1:length(ZR)
    for ji = 1:length(ZI)
        x = zr(jr); y = zi(ji);
        S = zeros(s,1);
        Shat = zeros(s+1,1);
        Shat(1) = 1i*y;
        for ell = 1:s
            S(ell) = x/(1-x*A(ell,ell))*...
                (1+[0,A(ell,1:ell-1)]*[0;S(1:ell-1)]+...
                Ahat(ell+1,1:ell)*Shat(1:ell));
            Shat(ell+1) = 1i*y/x*S(ell);
        end
        R(ji,jr) = 1+b'*S+bhat'*Shat; % growth factor
    end
end
contourf(zr,zi,1-abs(R),[0 0])
axis equal
