function absolute_stability_region_demo(methods)
% Plot selected regions, for example absolute_stability_region_demo(1:8).
% Methods 1--6 are Taylor, 7 is backward Euler, and 8 is trapezoid.

% Define a grid in the complex plane.
[zr, zi] = meshgrid(-5:0.1:2, -3:0.1:3);
z = zr + 1i * zi;

for method = methods
    % Compute the stability function.
    if method <= 6
        R = 0;
        for k = 0:method
            R = R + z.^k / factorial(k);
        end
        method_name = sprintf('Taylor method, order %d', method);
    elseif method == 7
        R = 1 ./ (1 - z);
        method_name = 'Backward Euler';
    else
        R = (1 + z / 2) ./ (1 - z / 2);
        method_name = 'Implicit trapezoid method';
    end

    % Plot the stable and unstable regions in a new figure.
    figure
    % Use finite contour levels: Octave cannot fill contours whose level
    % list contains infinities.  The binary mask gives the same stable /
    % unstable partition in both MATLAB and Octave.
    stable = double(abs(R) <= 1);
    contourf(zr, zi, stable, [-0.5, 0.5, 1.5], 'k')
    title(method_name)
    xlabel('Re(z)')
    ylabel('Im(z)')
    axis equal
    grid on
end
