% Script for the fully implicit method, using LU factorization.
% Initial values are to be provided by u0(x).
% Boundary values are to be provided by uL(x,tau), uH(x,tau).
%
% Before runing the script, m-files for u0, uL, uH need to
% be provided and values given to the following variables:
% xL = lower boundary
% xH = upper boundary
% N = number of spatial intervals (less 1)
% M = number of time steps
% T= final time
%
dx=(xH-xL)/(N+1);
dt=T/M;
alpha=dt/(dx^2);
e=ones(N,1);   %Note: this needs to be a column
A=spdiags([-alpha*e,(1+2*alpha)*e,-alpha*e],[-1,0,1],N,N);
[L,U]=lu(A);
bL=sparse(1,1,alpha,N,1);  %[alpha;zeros(N-1,1)];
bH=sparse(N,1,alpha,N,1);  %[zeros(N-1,1);alpha];
x=linspace(xL,xH,N+2);
u=u0(x(2:N+1))';
tau=0;
%
%Update Loop:
for m=1:M
     tau=m*dt;
     rhs=u+uL(xL,tau)*bL+uH(xH,tau)*bH;
     u=U\(L\rhs);
end
u=[uL(xL,T),u',uH(xH,T)];
