% heat_01.m
% Adapted from the code written initially by Dr. D. Pelinovsky
% (Assignment 5, Problem 3);


clear all;
close all;

% setting the grid in space;
h = 0.1; x = 0 : h : 1; n = length(x) - 1;

% setting the grid in time;
tau = 0.001; t = 0 : tau : 0.1; r = tau/h^2;
U1 = zeros(length(x),length(t));

% setting the initial condition and the RHS;
U1(:,1) = sin(pi*x');  U2 = U1; U3 = U1;
f = sin(3*pi*x)';

% assembling the matrices;
A1 = (1 - 2*r)*diag(ones(1,n-1)) + r*(diag(ones(1,n-2),1)+diag(ones(1,n-2),-1));
A2 = (1 + 2*r)*diag(ones(1,n-1)) - r*(diag(ones(1,n-2),1)+diag(ones(1,n-2),-1));
A3 = (1 - r)*diag(ones(1,n-1)) + 0.5*r*(diag(ones(1,n-2),1)+diag(ones(1,n-2),-1));
A4 = (1 + r)*diag(ones(1,n-1)) - 0.5*r*(diag(ones(1,n-2),1)+diag(ones(1,n-2),-1));

for k = 1 : (length(t)-1)
% the explicit method;  
    U1(2:n,k+1) = A1*U1(2:n,k) + tau*f(2:n);
% the implicit solution;    
    U2(2:n,k+1) = A2 \ (U2(2:n,k) + tau*f(2:n));
% the Crank-Nicolson solution;    
    U3(2:n,k+1) = A4 \ (A3*U3(2:n,k) + tau*f(2:n));
end

% the analytical solution;
Uex = exp(-pi^2*0.1)*sin(pi*x)+sin(3*pi*x)/(9*pi^2);


figure(1);
plot(x',U1(:,length(t)),':*r',x',U2(:,length(t)),':*b',x',U3(:,length(t)),':*g',x,Uex,'-c');
xlabel('x'); ylabel('u(x,t) at t = 0.1'); legend('explicit','implicit','CN');



% extending the integration interval to exhibit the instability
% of the explicit scheme;
tau
t = 0 : tau : 0.5; 
U = zeros(length(x),length(t));
U(:,1) = sin(pi*x');  
for k = 1 : (length(t)-1)
    U(2:n,k+1) = A1*U(2:n,k) + tau*f(2:n);
end
figure(2);
[T,X] = meshgrid(t,x);
surf(X,T,U); 
xlabel('x'); ylabel('t'); zlabel('u'); 



% plotting the CN solution profiles at various times;
t = 0 : tau : 1; 
U = zeros(length(x),length(t));
U(:,1) = sin(pi*x');  
for k = 1 : (length(t)-1)
    U(2:n,k+1) = A4 \ (A3*U(2:n,k) + tau*f(2:n));
end
figure(3);
hold on;
for k = 1 : 20 : 101
    plot(x,U(:,k));
end
xlabel('x'); ylabel('u'); hold off;

