close all; clear; clc;

%% Show the stability region of the Forward Euler method
h = 1;
t = (0:round(100./h)).*h; % Time axis
y0 = 1;

theta = linspace(0,pi,1e2);

rs = linspace(-3,3,32); is = linspace(-3,3,32);
H = zeros(numel(rs),numel(is)); % Preallocate
as = randperm(numel(H));

% Prepare figure
figure(1); clf; hold on;
ih = imagesc(rs,is,H);
plot(rs,zeros(size(rs)),'k--');
plot(zeros(size(is)),is,'k--');
plot(sin(theta)-1,cos(theta),'k','Linewidth',1);
plot(-sin(theta)-1,cos(theta),'k','Linewidth',1);
xlim([min(rs) max(rs)]); ylim([min(is) max(is)]);
axis square; colormap(flipud(plasma));
drawnow;

for ai = 1:numel(H)
    [ii,ri] = ind2sub(size(H),as(ai));
    lambda = (rs(ri)+1i*is(ii))./h;
    dydt = @(t,y) lambda*y;
    [~,exact] = ode23(dydt,t,y0);
    y = zeros(numel(t),1); y(1,:) = y0;
    for ti = 2:numel(t)
        dy = h.*dydt(t(ti),y(ti-1,:)); % Change in y is the current gradient
        y(ti,:) = y(ti-1,:)+dy; % Go forward in time
    end
    %plot(t,exact-y); drawnow;
    value = log10(rms(exact-y));
    H(ii,ri) = value;
    if ~mod(ai,round(numel(H)/20)); ih.CData = H; drawnow; end
end
ih.CData = H;