function [x,W] = Crank_Nicholson(Fun,alpha,l,m,T,N)
%
% Solving the heat equation with Crank Nicholson
%
% Fun: function at t = 0;
% 
[FunFcn,msg] = fcnchk(Fun,0);
if ~isempty(msg)
    error('InvalidFUN',msg);
end
k = T/N;
h = l/m;
lambda = (alpha^2)*k/(h^2);
x = (1:m-1)'*h; 
W = zeros(m-1,N+1);
W(:,1) = FunFcn(x);
for j = 2:N+1
    w = (1-lambda)*W(:,j-1) + (lambda/2)*[W(2:end,j-1);0] + (lambda/2)*[0;W(1:end-1,j-1)];
    e = ones(m-1,1);
    A = spdiags([-(lambda/2)*e, (1+lambda)*e,-(lambda/2)*e],-1:1, m-1,m-1);
    W(:,j) = A\w;
end


