%
%  lect3.m -- numerical instability example
%
%      y' = Q y    on t = [0,L]  
%
%      y(0) = A
%
%  explicit midpoint method
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  set parameters 

Q = -1.0;  A = 1.0;

N = 64;  L = 3;  dt = L/N;

t  = 0:dt:L;
ys = zeros(1,N+1);

%  initial condition & fwd euler step (change 0 to 1 for no osc'ns)

B = (Q*dt + sqrt(1 + 0*(Q*dt)^2)) * A;

ys(1) = A;  ys(2) = B;

%  iteration loop

for j=2:1:N
  ys(j+1) = ys(j-1) + 2*dt*Q * ys(j);
end 

y_ex = exp(Q*t);

error = ys-y_ex;  
disp(['final error = ',num2str(error(end))]);

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
%  solution plot

clf;subplot(2,1,1)
plot(t,ys,'k.')
hold on
plot(t,ys,'b')
%plot(t,y_ex,'r')
title('\bf exponential decay')
xlabel('\bf t-axis');  ylabel('\bf y-axis');

subplot(2,1,2)
plot(t,ys-y_ex,'k.')
hold on
plot(t,ys-y_ex,'b')
title('\bf error growth of computational mode')
%title('\bf error of pure decay mode (rigged)')
xlabel('\bf t-axis');  ylabel('\bf y_{comp}-y_{exact}');