%begin cranknicholson.m

clear;close;

% Set the number of grid points and build the cell-center
% grid.

N=input(' Enter N, cell number - ')
L=10;
h=L/N;

x=-.5*h:h:L+.5*h;
% Turn x into a column vector.
x=x';

%  Load the diffusion coefficient array (just 1 for now--
%  later you will load it with something more interesting.)
%  Make it a column vector.

D=ones(N+2,1);

% Load Dm with average values D(j-1/2) and Dp with D(j+1/2)
% Make them column vectors.

Dm=zeros(N+2,1);Dp=zeros(N+2,1);
Dm(2:N+1)=.5*(D(2:N+1)+D(1:N)); % average j and j-1
Dp(2:N+1)=.5*(D(2:N+1)+D(3:N+2)); % average j and j+1

% Initialize the temperature with a sine function.
T=sin(pi*x/L);

% Find the maximum of T for setting the plot frame.
Tmax=max(T);Tmin=min(T);

% Suggest a timestep by giving the time step for
% explicit stability, then ask for the time step.

fprintf(' Maximum explicit time step: %g \n',h^2/max(D))
tau = input(' Enter the time step - ')

tfinal=input(' Enter the total run time - ')

% Set the number of time steps to take.
nsteps=tfinal/tau;

% Define a useful constant.
coeff=tau/h^2 ;

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

% load A and B at interior points using
% colon commands

for j=2:N+1
   A(j,j-1)= -0.5*Dm(j)*coeff;
   A(j,j)  = 1+.5*(Dm(j)+Dp(j))*coeff;
   A(j,j+1)= -0.5*Dp(j)*coeff;

   B(j,j-1)=  0.5*Dm(j)*coeff;
   B(j,j)  = 1-.5*(Dm(j)+Dp(j))*coeff;
   B(j,j+1)=  0.5*Dp(j)*coeff;
end

% now load the boundary conditions for the
% case T(0)=0 and T(L)=0

A(1,1)=0.5;A(1,2)=0.5;B(1,1)=0.;
A(N+2,N+1)=0.5;A(N+2,N+2)=0.5;B(N+2,N+2)=0;


% This is the time advance loop.
for mtime=1:nsteps

t=mtime*tau; % define the time

r=B*T;
r(1)=0;r(N+2)=0; % set the right side vector r for T(0)=0 and T(L)=0
T=A\r; % do the linear solve

% Make a plot of the radial T(r) profile every once in a while.
   if(rem(mtime,5) == 0)
   plot(x,T)
   axis([0 L Tmin Tmax])
   pause(.1)
   end

end