clear all; clc;
format short g;

RK1 =@(z) 1 + z;
RK2 =@(z) 1 + z + (z.^2)/2;
RK3 =@(z) 1 + z + (z.^2)/2 + (z.^3)/6;
RK4 =@(z) 1 + z + (z.^2)/2 + (z.^3)/6 + (z.^4)/24;

BE1 =@(z) (1 - z).^-1;

% The rectangular window of the Argand plane where you want to see the contours
X_min = -03;  X_max = 03;
Y_min = -03;  Y_max = 03;
xgrids = linspace(X_min,X_max,1001);
ygrids = linspace(Y_min,Y_max,1001);

[X,Y] = ndgrid(xgrids, ygrids);
Z = complex(X,Y);
absRK1 = abs(RK1(Z));
absRK2 = abs(RK2(Z));
absRK3 = abs(RK3(Z));
absRK4 = abs(RK4(Z));
absBE1 = abs(BE1(Z));

fig1 = cfigure(20,20);
figure(fig1);
hold on;
grid on;
axis equal;
% set(gcf,'DefaultLineLineWidth',2);
[CH1,h1] = contour(X,Y,absRK1,[1,1],'LineColor','k','LineWidth',1);
[CH2,h2] = contour(X,Y,absRK2,[1,1],'LineColor','b','LineWidth',1);
[CH3,h3] = contour(X,Y,absRK3,[1,1],'LineColor','r','LineWidth',1);
[CH4,h4] = contour(X,Y,absRK4,[1,1],'LineColor','m','LineWidth',1);
% [CH5,h5] = contour(X,Y,absBE1,[1,1],'LineColor','m');

% Get the normalized x and y coordinates and pack it using deal.
pos = get(gca, 'Position');
nx =@(aB,aE) deal(pos(1) + pos(3)*([aE(1) aB(1)]+ abs(min(xlim)))/diff(xlim));
ny =@(aB,aE) deal(pos(2) + pos(4)*([aE(2) aB(2)]+ abs(min(ylim)))/diff(ylim));

aB = [-.1,2.3]; aE = [0.5,2.3]; aStr = 'stability region of RK4';
annotation('textarrow',nx(aB,aE),ny(aB,aE),'String',aStr,'FontSize',14);

aB = [-.2,1.7]; aE = [0.5,1.7]; aStr = 'stability region of RK3';
annotation('textarrow',nx(aB,aE),ny(aB,aE),'String',aStr,'FontSize',14);

aB = [-.4,1.1]; aE = [0.5,1.1]; aStr = 'stability region of RK2';
annotation('textarrow',nx(aB,aE),ny(aB,aE),'String',aStr,'FontSize',14);

aB = [-.5,0.5]; aE = [0.5,0.5]; aStr = 'stability region of RK1';
annotation('textarrow',nx(aB,aE),ny(aB,aE),'String',aStr,'FontSize',14);

view(0,90);
set(gca,'XLim',[X_min X_max]);
set(gca,'YLim',[Y_min Y_max]);
set(gca,'XTick',X_min:1:X_max);
set(gca,'YTick',Y_min:1:Y_max);
ylabel('Im(\lambda h)','FontSize',16);
xlabel('Re(\lambda h)','FontSize',16);
title_str = ['Absolute stability regions of some RK methods for u''=\lambda u'];
title(title_str,'FontSize',16);
hold off;


