% homogeneous Dirichlet boundary conditions.  First, I compute with a
% time-step that is small enough for Forward Euler to be numerically
% stable.  I plot the solution at various times, I plot the log of the norm
% of the solution at various times, and I plot the power spectrum at
% various times.
%
% I then repeat this with a larger time step for which Forward Euler is
% numerically unstable.
%
% I then repeat this, with the same time step, using Crank-Nicolson for
% which there is no stability constratint on the time step.


t_0 = 0;
t_f = .5;
N = 128;
M = 1000;

[u1,err,x,t1] = heat_eul_dir(t_0,t_f,M,N);
figure(1)
for j=0:10
    plot(x,u1(:,j*100+1),'.-')
    title(t1(100*j+1),'FontSize',16)
    xlabel('Explicit Euler','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))))
title('log10(maximum of u) versus time','FontSize',16)


display('hit any key for more')
pause

figure(1)
for j=0:10
    [k,amp] = find_spec(u1(1:N,100*j+1));
    plot(k,amp,'.-')
    title(t1(100*j+1),'FontSize',16)
    axis([0,N/2,-17,2])
    pause(1)
end

display('hit any key for more')
pause

M = 800;
[u2,err,x,t2] = heat_eul_dir(t_0,t_f,M,N);
figure(1)
for j=0:10
    plot(x,u2(:,80*j+1),'.-')
    title(t2(80*j+1),'FontSize',16)
    xlabel('Forward Euler','FontSize',16)
    axis([0,2*pi,-1.1,1.1])
    pause(1)
end

figure(2)
plot(t2,log10(max(abs(u2))))
title('log10(maximum of u) versus time','FontSize',16)
xlabel('notice that things start to grow exponentially shortly before time 0.3')

display('hit any key for more')
pause

figure(1)
for j=0:80
    [k,amp] = find_spec(u2(1:N,10*j+1));
    plot(k,amp,'.-')
    title(t2(10*j+1),'FontSize',16)
    axis([0,N/2,-17,2])
    pause(1)
end

display('hit any key for more')
pause

[u3,err,x,t1] = heat_cn(t_0,t_f,M,N);
figure(1)
for j=0:20
    plot(x,u3(:,j+1),'.-')
    axis([0,2*pi,-1.1,1.1])
    title(t1(j+1),'FontSize',16)
    xlabel('Crank-Nicolson','FontSize',16)
    pause(1)
end

figure(2)
plot(t1,log10(max(abs(u3))))
title('log10(maximum of u) versus time','FontSize',16)

display('hit any key for more')
pause

figure(1)
for j=0:20
    [k,amp] = find_spec(u3(1:N,j+1));
    plot(k,amp,'.-')
    axis([0,N/2,-17,2])
    pause(1)
end


