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 = 1001;
h  = L/Nx;
x  = linspace(0,L,Nx+1);

% Set timestep size 
dt = 0.001;

% Calculate number of iterations (Nmax) needed to iterate to t=Tmax
Tmax = 0.5;
Nmax = floor(Tmax/dt);
t    = linspace(0,Nmax*dt,Nmax+1);
Tmid = zeros(Nmax,1);
Tmid(1) = Tinit;

% 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);

figure(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;

  % 2nd order bd
%  if (n==1),
%    v1 = v;
%    v = (I - dt*A)\(v1 + dt*b);
%  else,
%    v0 = v1;
%    v1 = v;
%    v = (I - 2/3*dt*A)\(4/3*v1 - 1/3*v0 + 2/3*dt*b);
%  end
    
  T = [Tleft; v; Tright];
  plot(x,T,'*'); 
  axis([0,1,400,1200]);
  drawnow;

  Tmid(n+1) = T(floor(0.5*(Nx+1)));
end


figure(2);
plot(t,Tmid);hold on;

figure(3);
plot(x,T); hold on;
axis([0,1,400,1200]);