% This code solve the one-dimensional heat diffusion equation
% for the problem of a bar which is initially at T=Tinit and
% suddenly the temperatures at the left and right change to
% Tleft and Tright.
%
% Upon discretization in space by a finite difference method,
% the result is a system of ODE's of the form,
%
% u_t = Au + b
%
% The code calculates A and b.  Then, these are used in
% various different time integration methods to evolve the
% temperature in time.  See lecture notes for more info.
%
% Using this code, we will talk about the following issues:
%
% * Use of explicit methods for stiff systems
% * Implementation of implicit methods for linear systems
% * Sparse vs. full matrices
% * Behavior of trapezoidal when lambda dt -> -infinity.

clear all;
close all;

% Set non-dimensional thermal coefficient
k   = 1.0; % this is really k/(rho*cp)

% Set length of bar
L = 1.0; % non-dimensional

% Set initial temperature
Tinit = 400;

% Set left and right temperatures for t>0
Tleft = 800;
Tright = 1000;

% Set up grid size
Nx = 1000;
h  = L/Nx;
x  = linspace(0,L,Nx+1);

% Set timestep size 
dt = 1e-1;

% Calculate number of iterations (Nmax) needed to iterate to t=Tmax
Tmax = 0.5;
Nmax = floor(Tmax/dt);

% Initialize a sparse matrix to hold stiffness & identity matrix
A = spalloc(Nx-1,Nx-1,3*(Nx-1));
I = speye(Nx-1);

% These commands allocate full (non-sparse) matrices
%A = zeros(Nx-1,Nx-1);
%I = eye(Nx-1);

% Calculate stiffness matrix

for ii = 1:Nx-1,
    
  if (ii > 1),
    A(ii,ii-1) = k/h^2;
  end
  
  if (ii < Nx-1),
    A(ii,ii+1) = k/h^2;
  end
  
  A(ii,ii) = -2*k/h^2;
  
end

% Set forcing vector
b = zeros(Nx-1,1);
b(1)    = k*Tleft/h^2;
b(Nx-1) = k*Tright/h^2;

% Set initial vector
v = Tinit*ones(Nx-1,1);

% Iterate in time 
for n = 1:Nmax,
  
  % Backward Euler
  %v = (I - dt*A)\(v + dt*b);

  % Forward Euler
  %v = v + dt*(A*v+b);

  % Trapezoidal
  rhs = v + 0.5*dt*A*v + dt*b;
  G = I - 0.5*dt*A;
  v = G\rhs;

  T = [Tleft; v; Tright];
  plot(x,T,'*'); 
  axis([0,1,400,1200]);
  drawnow;
  
end