function temple5044_absolute_stability_region_rk
%TEMPLE5044_ABSOLUTE_STABILITY_REGION_RK
%   Plots the regions of absolute stability for a number
%   of Runge-Kutta schemes given by their Butcher tableau.
%
%   (C) 2019/03/15 by Benjamin Seibold
%       http://www.math.temple.edu/~seibold/

%========================================================================
% Define RK schemes
%========================================================================
name{1} = 'forward Euler';
A{1} = [0];
b{1} = [1];

name{end+1} = 'Heun''s method';
A{end+1} = [0 0;1 0];
b{end+1} = [1/2;1/2];

name{end+1} = 'Shu-Osher method';
A{end+1} = [0 0 0;1 0 0;1/4 1/4 0];
b{end+1} = [1/6;1/6;2/3];

name{end+1} = 'RK4';
A{end+1} = [0 0 0 0;1/2 0 0 0;0 1/2 0 0;0 0 1 0];
b{end+1} = [1;2;2;1]/6;

name{end+1} = 'backward Euler';
A{end+1} = [1];
b{end+1} = [1];

name{end+1} = 'Crank-Nicolson';
A{end+1} = [0 0;1/2 1/2];
b{end+1} = [1/2;1/2];

name{end+1} = 'Gauss-Legendre 4';
A{end+1} = [1/4 1/4-sqrt(3)/6;1/4+sqrt(3)/6 1/4];
b{end+1} = [1/2;1/2];

name{end+1} = 'TR-BDF2';
A{end+1} = [0 0 0;1/4 1/4 0;1/3 1/3 1/3];
b{end+1} = [1/3;1/3;1/3];

name{end+1} = 'Radau IA 2';
A{end+1} = [1/4 -1/4;1/4 5/12];
b{end+1} = [1/4;3/4];

%========================================================================
% Plot stability regions
%========================================================================
ax = 8; % 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},ax) % plot RAS for that scheme
    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,ax)
if nargin<3, ax = 1; end % size of stability region plot
s = length(b); % number of stages
e = ones(s,1); % vector of all ones
I = eye(s); % identity matrix
zr = linspace(-1,1,100)*ax; zi = zr;
[ZR,ZI] = meshgrid(zr,zi);
R = ZR*0;
for jr = 1:length(ZR)
    for ji = 1:length(ZI)
        z = zr(jr)+1i*zi(ji);
        R(ji,jr) = 1+z*b'*((I-z*A)\e);
    end
end
contourf(zr,zi,1-abs(R),[0 0])
axis equal
