% Plot the stabilty domain for a given stability function. 
% Domain: [-a, a, -b, b]
a = 4
b = 4
[x, y] = meshgrid(linspace(-a, a), linspace(-b, b));
z = x + i * y;

% Stability function.
R = abs(1 + z + z.^2 / 2 + z.^3 / 6);                                       % RK3
%R = abs(1 + z + z.^2 / 2 + z.^3 / 6 + z.^4 / 24 + z.^5 / 120 + z.^6 / 600); % DOPRI5

% Make the plot.
contourf(x, y, R, [1, 1], 'k')
axis equal
axis([-a, a, -b, b])
grid on
hold on
plot([-a, a], [ 0, 0], 'k', 'LineWidth', 1);
plot([ 0, 0], [-a, a], 'k', 'LineWidth', 1);
