% stabilityregion Plot the stability region of Runge-Kutta methods up to
%                 order 4
% Lei Zhang, 25-09-2012
function stabilityregion(choice)

clear i;
switch choice
    case 1 % whole region
        [X,Y] = meshgrid(-5:0.005:5,-5:0.005:5);
    case 2 % close to 0
        [X,Y] = meshgrid(-0.05:5e-5:0.05,-0.05:5e-5:0.05);
end
Mu = X+i*Y;
R1 = 1 + Mu;
R2 = R1 + .5*Mu.^2;
R3 = R2 + (1/6)*Mu.^3;
R4 = R3 + (1/24)*Mu.^4;
Rhat1 = abs(R1);
Rhat2 = abs(R2);
Rhat3 = abs(R3);
Rhat4 = abs(R4);
figure; axis equal; width = 2;
switch choice
    case 1
        axis([-4 4 -4 4]);
    case 2
        axis([-0.05 0.05 -0.05 0.05]);
end
hold on
[C,h]=contour(X,Y,Rhat1,[1, 1]); % plot the contour line for abs(R) = 1
set(h,'Color','b','linewidth',width)
pause
[C,h]=contour(X,Y,Rhat2,[1, 1]);
set(h,'Color','g','linewidth',width)
pause
[C,h]=contour(X,Y,Rhat3,[1, 1]);
set(h,'Color','r','linewidth',width)
pause
[C,h]=contour(X,Y,Rhat4,[1, 1]);
set(h,'Color','c','linewidth',width)
legend('RK1','RK2','RK3','RK4')
line([-4 4],[0 0]);
line([0 0],[-4 4]);
title('Stability regions of Runge-Kutta methods up to order 4')