  function uNext = CrankN(uNow,a,b,n,tc,deltaT,s)
% uNow is a column n-vector.
% a < b 
% n is a positive integer
% tc is the ``current time''.
% deltaT >0 is the time step. 
% s is a handle to a function of the form s(x,t)
% uNext is a column n-vector whose entries satisfy uNext(1) = uNow(1),
%     uNext(n) = uNow(n), and 


%  uNew(k) - uNow(k)      1   (uNew(k+1)-2*uNew(k)+uNew(k-1)) 
%  -----------------  =  ---  -----------------------------   + 
%      deltaT             2                h^2

%                         1   (uNow(k+1)-2*uNow(k)+uNow(k-1)) 
%                        ---  -----------------------------   + s(x(k),tc+deltaT/2)
%                         2                h^2

% for k=2:n-1 where 
h = (b-a)/(n-1);
x = linspace(a,b,n);


r = deltaT/(2*h^2);
% then we have a set of n-2 linear equations:

%   -r*uNew(k-1) + (1+2r)*uNew(k) - r*uNew(k+1)  = deltaT*s(x(k),tc+deltaT/2) 
%    r*uCur(k-1) + (1-2r)*uCur(k) - r*uCur(k+1)

 
 uNext = zeros(n,1);
 uNext(1) = uNow(1);
 uNext(n) = uNow(n);
 % Set up a system for uNext(2:n-1)...
 T = zeros(n-2,n-2);
 rhs = zeros(n-2,1);
 for k=1:n-2
     T(k,k) = 1+2*r;
     if k>1
         T(k,k-1) = -r;
     end
     if k<n-2
         T(k,k+1) = -r;
     end
     rhs(k) = deltaT*s(x(k+1),tc+deltaT/2) + r*uNow(k+2) + (1-2*r)*uNow(k+1) + r*uNow(k);
 end
 rhs(1) = rhs(1)+r*uNext(1);
 rhs(n-2)=rhs(n-2)+r*uNext(n);
 uNext(2:n-1) = T\rhs;