% 
% Specify range of S= sigma*dt

smax = 2.5;
ds = 0.1;
sr = -smax:ds:smax;
si = -smax:ds:smax;
[Sr Si] = meshgrid(sr,si);
S = Sr + Si*1i;

% Set up red-blue color map

cred = [1 0 0];
cwhite= [1 1 1];
cblue = [0 0 1];
cfrac = (0.1:0.1:1)';
cbwmap = flipud(cfrac)*cblue + flipud(1-cfrac)*cwhite;
cwrmap = cfrac*cred + (1-cfrac)*cwhite;
cbrmap = [cbwmap;cwhite;cwrmap];

subplot(2,2,1)
A = exp(S);

pcolor(sr,si,log(abs(A)))
caxis([-5 5])
colormap(cbrmap)
shading('interp')
text(-1,2,'Color:log|A|')
colorbar
hold on

contour(sr,si,abs(A),[1,1],'k');
xlabel('Re[\sigma\Delta t]')
ylabel('Im[\sigma\Delta t]')
text(0.5,0,'|A|>1')
text(-1.5,0,'|A|<1')
plot([0 0],[-2.5 2.5],'k--')
axis square
title('Exact')

subplot(2,2,2)
A = 1 + S;

pcolor(sr,si,log(abs(A)))
caxis([-5 5])
colormap(cbrmap)
shading('interp')
text(-1,2,'Color:log|A|')
colorbar
hold on

contour(sr,si,abs(A),[1,1],'k');
xlabel('Re[\sigma\Delta t]')
ylabel('Im[\sigma\Delta t]')
text(0.5,0,'|A|>1')
text(-1.5,0,'|A|<1')
plot([0 0],[-2.5 2.5],'k--')
axis square
title('Forward Euler')

subplot(2,2,3)
A = 1./(1 - S);

pcolor(sr,si,log(abs(A)))
caxis([-5 5])
colormap(cbrmap)
shading('interp')
text(-1,2,'Color:log|A|')
colorbar
hold on

contour(sr,si,abs(A),[1,1]);
xlabel('Re[\sigma\Delta t]')
ylabel('Im[\sigma]\Delta t]')
text(0.5,0,'|A|>1')
text(-1.5,0,'|A|<1')
plot([0 0],[-2.5 2.5],'k--')
axis square
title('Backward Euler')

subplot(2,2,4)
A = (1 + S/2)./(1 - S/2);

pcolor(sr,si,log(abs(A)))
caxis([-5 5])
colormap(cbrmap)
shading('interp')
text(-1,2,'Color:log|A|')
colorbar
hold on

contour(sr,si,abs(A),[1,1],'k');
xlabel('Re[\sigma\Delta t]')
ylabel('Im[\sigma\Delta t]')
text(0.5,0,'|A|>1')
text(-1.5,0,'|A|<1')
plot([0 0],[-2.5 2.5],'k--')
axis square
title('Trapezoidal')
