%
% Plots absolute-stability regions for explicit RK methods of orders 1-8. 
% 

p = 8;

xmin = -5; xmax = 3; 
ymin = -5; ymax = 5;

hx = 100\(xmax - xmin);
hy = 100\(ymax - ymin);

x = xmin:hx:xmax; y = ymin:hy:ymax; 
[X,Y]=meshgrid(x,y);
Z = complex(X,Y); 

figure(1); clf; hold on; 
axis([xmin xmax ymin ymax]);grid on;axis equal;

plot([xmin xmax],[0 0],'-k','LineWidth',1);
plot([0 0],[ymin ymax],'-k','LineWidth',1);

ASRT = 'Absolute Stability Regions for Explicit Runge-Kutta Methods';
title(ASRT,'FontSize',16);

% Use MATLAB's ColorOrder colors.
clrs = [0    0.4470    0.7410;
    0.8500    0.3250    0.0980;
    0.9290    0.6940    0.1250;
    0.4940    0.1840    0.5560;
    0.4660    0.6740    0.1880;
    0.3010    0.7450    0.9330;
    0.6350    0.0780    0.1840];

for k = 1:p
    switch k
        case 1
            range = [-4 2 -4 4];
            g = @(Z)1+Z;
            ABSGVALS= abs(feval(g,Z));
            contour(X,Y,ABSGVALS,[1 1+10*eps],'LineWidth',2,...
                'LineColor',clrs(1,:) );
            title({ASRT;'Order 1 (Forward Euler)'});
            pause;
        case 2
            range = [-4 2 -4 4];
            g = @(Z)1+Z+2\Z.^2;
            ABSGVALS= abs(feval(g,Z));
            contour(X,Y,ABSGVALS,[1 1+10*eps],'LineWidth',2,...
                'LineColor',clrs(2,:) );
            title({ASRT;'Orders 1-2'});
            pause;
        case 3 
            range = [-4 2 -4 4];
            g = @(Z)1+Z+2\Z.^2+6\Z.^3;
            ABSGVALS= abs(feval(g,Z));
            contour(X,Y,ABSGVALS,[1 1+10*eps],'LineWidth',2,...
                'LineColor',clrs(3,:) );
            title({ASRT;'Orders 1-3'});
            pause;
        case 4
            range = [-4 2 -4 4];
            g = @(Z)1+Z+2\Z.^2+6\Z.^3+24\Z.^4;
            ABSGVALS= abs(feval(g,Z));
            contour(X,Y,ABSGVALS,[1 1+10*eps],'LineWidth',2,...
                'LineColor',clrs(4,:) );
            title({ASRT;'Orders 1-4'});
            pause;
        case 5
            range = [-5 2 -5 5];
            g = @(Z)1+Z+2\Z.^2+6\Z.^3+24\Z.^4+120\Z.^5;
            ABSGVALS= abs(feval(g,Z));
            contour(X,Y,ABSGVALS,[1 1+10*eps],'LineWidth',2,...
                'LineColor',clrs(5,:) );
            title({ASRT;'Orders 1-5'},'FontSize',16);
            pause;
        case 6
            range = [-5 2 -5 5];
            g = @(Z)1+Z+2\Z.^2+6\Z.^3+24\Z.^4+120\Z.^5+720\Z.^6;
            ABSGVALS= abs(feval(g,Z));
            contour(X,Y,ABSGVALS,[1 1+10*eps],'LineWidth',2,...
                'LineColor',clrs(6,:) );
            title({ASRT;'Orders 1-6'});
            pause;
        case 7
            range = [-5 2 -5 5];
            g = @(Z)1+Z+2\Z.^2+6\Z.^3+24\Z.^4+120\Z.^5+720\Z.^6+5040\Z.^7;
            ABSGVALS= abs(feval(g,Z));
            contour(X,Y,ABSGVALS,[1 1+10*eps],'LineWidth',2,...
                'LineColor',clrs(7,:) );
            title({ASRT;'Orders 1-7'});
            pause;
        case 8
            range = [-5 3 -5 5];
            g = @(Z)1+Z+2\Z.^2+6\Z.^3+24\Z.^4+120\Z.^5+720\Z.^6+...
                5040\Z.^7+40320\Z.^8;
            ABSGVALS= abs(feval(g,Z));
            contour(X,Y,ABSGVALS,[1 1+10*eps],'LineWidth',2,...
                'LineColor',clrs(1,:) );
            title({ASRT;'Orders 1-8'});
            
    end
              
end

hold off;