% This is a demo of Forward Euler time-stepping for u_t = D u_xx with both
% homogeneous Dirichlet and homogeneous Neumann boundary conditions.  
%
% First, I consider the homogeneous Dirichlet boundary conditions.  I take
% the time-step is "too large" (lambda = 8.3) and plot the true solution 
% along with the discrete approximation at various timesteps.  I see
% instability, as expected.  Also, at each time I compute the L^\infty
% norm of the discrete solution and plot the log of this versus time.  For
% the true solution this would lead to a linear plot with constant negative
% slope.  For the unstable computation one sees this at first but then the
% growth due to round-off error appears.
%
% I then repeat all this with a "small enough" time step (lambda = 0.42).  
%
% Finally, I redo both cases for the homogeneous Neumann boundary
% conditions.
help heat_eul_dir

display('hit any key for more')
pause

t_0 = 0;
t_f = 1;
N = 128;
M = 100;

[u1,err1,x,t1] = heat_eul_dir(t_0,t_f,M,N);
figure(1)
for j=0:20
     plot(x,u1(:,j+1),'o-')
     hold on
     plot(x,u1(:,j+1)-err1(:,j+1),'r')
     figure(1)
     hold off
     title(t1(j+1),'FontSize',16)
     xlabel('Explicit Euler, circles are approx. solution','FontSize',16)
     axis([0,2*pi,-1.1,1.1])
     pause(1)
end

display('hit any key for more')
pause

figure(2)
plot(t1,log10(max(abs(u1))),'o-')
axis([0,.25,-1,5])
title('log10(maximum of u) versus time','FontSize',16)
xlabel('notice that things start to grow exponentially shortly before time 0.2')

display('hit any key for more')
pause


M = 2000;
[u2,err2,x,t2] = heat_eul_dir(t_0,t_f,M,N);
figure(2)
for j=0:20
     plot(x,u2(:,10*j+1),'o-')
     hold on
     plot(x,u2(:,10*j+1)-err2(:,10*j+1),'r')
     figure(2)
     hold off
     axis([0,2*pi,-1.1,1.1])
     xlabel('Explicit Euler, circles are approx. solution','FontSize',16)
     title(t2(10*j+1),'FontSize',16)
     pause(1)
end

figure(2)
plot(t2,log10(max(abs(u2))),'o-')
axis([0,.25,-1,0.1])
title('log10(maximum of u) versus time','FontSize',16)

display('hit any key for more')
pause


help heat_eul_neu

display('hit any key for more')
pause


M = 100;
[u3,err3,x,t3] = heat_eul_neu(t_0,t_f,M,N);
figure(1)
for j=0:20
     plot(x,u3(:,j+1),'o-')
     hold on
     plot(x,u3(:,j+1)-err3(:,j+1),'r')
     hold off
     figure(1)
     axis([0,2*pi,-1.1,1.1])
     xlabel('fully implicit, circles are approx. solution','FontSize',16)
     title(t1(j+1),'FontSize',16)
     pause(1)
end

display('hit any key for more')
pause

figure(1)
plot(t3,log10(max(abs(u3))),'o-')
axis([0,.25,-1,5])
title('log10(maximum of u) versus time','FontSize',16)
xlabel('notice that things start to grow exponentially shortly before time 0.2')

display('hit any key for more')
pause


M = 2000;
[u4,err4,x,t4] = heat_eul_dir(t_0,t_f,M,N);
figure(2)
for j=0:20
     plot(x,u4(:,10*j+1),'o-')
     hold on
     plot(x,u4(:,10*j+1)-err4(:,10*j+1),'r')
     figure(2)
     hold off
     axis([0,2*pi,-1.1,1.1])
     xlabel('Explicit Euler, circles are approx. solution','FontSize',16)
     title(t2(10*j+1),'FontSize',16)
     pause(1)
end

display('hit any key for more')
pause

figure(2)
plot(t4,log10(max(abs(u4))),'o-')
axis([0,.25,-.7,0.1])
title('log10(maximum of u) versus time','FontSize',16)

