% Q2.2 Module 4
%
% Code to solve the 1D heat equation implicitly using Crank-Nicolson
% Student name: Renee Meaney
% Student ID: 27848779


clear all;

% Set the number of grid points
M=100; % These are the interior points (There will be M+2 including B.C's)
L=1.0;
dx=L/M;
x=-0.5*dx:dx:L+0.5*dx; % x-grid including outer BC points
x=x';  % Turn x into a column vector.
a=0.4;

% Setup the time domain.
dt = 0.01; % Timestep
tfinal = 1;
N = tfinal/dt;
t = linspace(0,tfinal,N+1);

% Introduce the initial condition via f(x)
f = zeros(M,1);
for m = 1:M
    if x(m) >= 0 & x(m) < L/4.0
        f(m) = 0; end
    if x(m) >= L/4 & x(m) < 3*L/4.0
        f(m) = sin((2*pi/L)*(x(m)-(L/4.0))); end   % Basic sinusoidal I.C.
    if x(m) >= 3*L/4.0 & x(m) <= L
        f(m) = 0; end
end
u = zeros(M+2,N);
u(:,1)=[0;f;0];
% So u(:,2) will be the solution at t = dt, u(:,3) at t = 2*dt etc.

% Create the matrices A and B by loading them with zeros
A=zeros(M+2);
B=zeros(M+2);

% load A and B at interior points
alpha = (1/2)*(a^2/dx^2);
beta =  (1/dt)+(a^2/dx^2);
gamma = (1/dt)-(a^2/dx^2);
for j=2:M+1
  A(j,j-1)= -alpha;
  A(j,j)  = beta;
  A(j,j+1)= -alpha;
  B(j,j-1)= alpha;
  B(j,j)  = gamma;
  B(j,j+1)= alpha;
end

% Add values of '1' to A(1,1) and A(M+2,M+2)
A(1,1)=1.0;            % T(0)=0
A(M+2,M+2)=1.0;        % T(L)=0

% Now we advance through time, solving for each new solution vector!
for n=1:N

u(:,n+1)= (A)\(B*u(:,n)); %Calculate the new solution vector u

% Make a movie of u in the loop, like in the previous module (with labels!)
% You should delete this before submitting and only generate Fig1 later.
%figure(1); %-%
%plot(x,u(:,n+1)) %-%
%xlabel('x') %-%
%ylabel('u(x,t)') %-%
%axis([0 1 0 1]) %-%
%pause(0.02); %-%

end


Fig1 = figure;


% Make a surface plot of u(x,t) using the 'surf' or 'mesh' functions.
% You do remember how to use 'meshgrid' right?

[x,t] = meshgrid(x,t);
surf(x,t,u');
title('Heat Equation surface plot using the Crank-Nicholson scheme');
xlabel('x');
ylabel('t');
zlabel('u(x,t)');

%Conclusions: