% Homework 12 - Bonus Problem:
%  Solving the Heat equation with zero boundary conditions using the
%  modified Euler's method.

clear all

r=0.6;
% K=0.006;      %  step size in t
h=0.1;        %  step size in  x
K = r*h^2;
% r=K/h^2; 
tmax=0.5;   

x=0:h:1;
M=length(x)-1;  
t=0:K:tmax;
N=length(t);    % THIS N = [(N+1) OF THE NOTES], i.e. N_notes = thisN - 1.

% Set up the boundary conditions for both U and the auxiliary variable Ubar:
U(1,:)=zeros(1,N);
U(M+1,:)=zeros(1,N);
U(:,1)=sin(pi*x');
Ubar(1)=0;
Ubar(M+1)=0;

% Advance the solution from the n-th to the (n+1)-st level.
%  First, compute the Ubar FOR ALL m,  and then compute U(m,n+1).
for n=1:N-1
    for m=2:M
        Ubar(m)=r*U(m+1,n)+(1-2*r)*U(m,n)+r*U(m-1,n);  % this is eq. (12.12)
        %  V V V V V V 
        % Surprisingly, if one computes Ubar using (12.17) instead of (12.12),
        %  it changes the round-off error! (As long as the round-off error
        %  does not become visible in the numerical error, the results of 
        %  using either 12.12 or 12.17 for computing Ubar are the same.)
        % ^ ^ ^ ^ ^ ^ ^
    end
    for m=2:M
        U(m,n+1) = U(m,n) + r/2*( (U(m+1,n)-2*U(m,n)+U(m-1,n)) + ...
                   (Ubar(m+1)-2*Ubar(m)+Ubar(m-1)) );
    end
    figure(120501);
    plot(x,sin(pi*x)*exp(-pi^2*t(n+1)),'r:',x,U(:,n+1))
    legend('exact','approx');
    xlabel('X', 'fontsize', 12)
    title(['Solution at t = ', num2str(t(n+1))], 'fontsize', 13)
    pause(0.2)
end

% Exact solution:
uexact=sin(pi*x)*exp(-pi^2*t(end));

figure(120502);
subplot(2,1,1);
plot(x,uexact,':',x,U(:,end));
legend('exact','approx');
xlabel('X', 'fontsize', 13)
ylabel('Solution at T_{max}', 'fontsize', 13)
subplot(2,1,2);
plot(x,uexact-U(:,end)');
xlabel('X', 'fontsize', 13)
ylabel('Error at T_{max}', 'fontsize', 13)
