% crank-nicholson scheme for 1d diffusion eqn. with fixed BC
function M = cn;

D = 4;
dt = 1;
dx = 1;
time = 1000;
length = 10;
s = 0; 

Ncells = length / dx + 1; 
Nsteps = time / dt + 1;
L = D*dt/((dx)^2)
  
u(Nsteps, Ncells)=0;
u(1,:) = 25;       % initial condition (time=0)
u(:,1) = 0;        % left BC  (x=0)
u(:,Ncells) = 100; % right BC (x=L)

diag_M = (2*(1-s)*L+1)*ones(Ncells,1);
diag_M(1) = 1;
diag_M(Ncells) = 1;
M = diag( diag_M );
M = M + diag( (s-1)*L*ones(Ncells-1,1),-1) + diag( (s-1)*L*ones(Ncells-1,1),1);
M( 1, 2 ) = 0;
M( Ncells, Ncells-1 ) = 0;
inv_M = inv( M );

for i = 1:Nsteps
  b(1)=0;
  b(Ncells)=100;
  for j = 2:(Ncells-1)
    b(j) = s*L*(u(i,j-1)+u(i,j+1)) + (1-2*s*L)*u(i,j);
  end
  u(i+1,:) = ( inv_M * b' )';
end

%figure;
%surf(u(1:(Nsteps/12.5):Nsteps,1:(Ncells/12.5):Ncells));
%view(-154,26);
figure;
plot(u(1:(Nsteps-1)/40:Nsteps,:)');