function [x,t,u] = MOL_BEuler(d0,d1,f,G,T,h,m)

% Use the method of lines to solve
%    u_t = u_xx + G(x,t), 0 < x < 1, 0 < t < T
% with boundary conditions
%    u(0,t) = d0(t), u(1,t) = d1(t)
% and initial condition
%    u(x,0) = f(x).
% Use the backward Euler's method to solve the system of
% ODEs. For the discretization, use a spatial stepsize of
% delta=1/m and a time step of h. 
%

x = linspace(0,1,m+1)'; delta = 1/m; delta_sqr = delta^2;
t = (0:h:T)'; N = length(t);

% Initialize u.
u = zeros(m+1,N);
u(:,1) = f(x);
u(1,:) = d0(t); u(m+1,:) = d1(t);
% Create tridiagonal coefficient matrix.
a = -(h/delta_sqr)*ones(m-1,1); c = a; 
b = (1+2*h/delta_sqr)*ones(m-1,1);
a(1) = 0; c(m-1) = 0; option = 0;

% Solve for u using the backward Euler's method.
for n=2:N
    g = G(x(2:m),t(n));
    g(1) = g(1) + (1/delta_sqr)*u(1,n);
    g(m-1) = g(m-1) + (1/delta_sqr)*u(m+1,n);
    f = u(2:m,n-1) + h*g;
    switch option
        case 0
            [v,alpha,beta,message] = tridiag(a,b,c,f,m-1,option);
            option = 1;
        case 1
            v = tridiag(alpha,beta,c,f,m-1,option);
    end
    u(2:m,n) = v;
end
u = u';
end % MOL_BEuler

function [x, alpha, beta, message] = tridiag(a,b,c,f,n,option)

% function [x, alpha, beta, message] = tridiag(a,b,c,f,n,option)
%
% Solve a tridiagonal linear system M*x=f
%
% INPUT:
% The order of the linear system is given as n.
% The subdiagonal, diagonal, and superdiagonal of M are given
% by the arrays a,b,c, respectively.  More precisely,
%     M(i,i-1) = a(i), i=2,...,n
%     M(i,i)   = b(i), i=1,...,n
%     M(i,i+1) = c(i), i=1,...,n-1
% option=0 means that the original matrix M is given as 
%   specified above.
% option=1 means that the LU factorization of M is already
%   known and is stored in a,b,c.  This will have been 
%   accomplished by a previous call to this routine.  In that 
%   case, the vectors alpha and beta should have been 
%   substituted for a and b in the calling sequence.
% All input values are unchanged on exit from the routine.
%
% OUTPUT:
% Upon exit, the LU factorization of M is already known and
% is stored in alpha,beta,c.  The solution x is given as well.
% message=0 means the program was completed satisfactorily.
% message=1 means that a zero pivot element was encountered 
%   and thesolution process was abandoned. This case happens 
%   only when option=0.

if option == 0
    alpha = a; beta = b;
    alpha(1) = 0;
   % Compute LU factorization of matrix M.
   for j=2:n
      if beta(j-1) == 0
         message = 1; return
      end
      alpha(j) = alpha(j)/beta(j-1); 
      beta(j) = beta(j) - alpha(j)*c(j-1);
   end
   if beta(n) == 0 
      message = 1; return
   end
end

% Compute solution x to M*x = f using LU factorization of M.
% Do forward substitution to solve lower triangular system.
if option == 1
    alpha = a; beta = b;
end
x = f; message = 0;
for j=2:n
   x(j) = x(j) - alpha(j)*x(j-1);
end

% Do backward substitution to solve upper triangular system.
x(n) = x(n)/beta(n);
for j=n-1:-1:1
   x(j) = (x(j) - c(j)*x(j+1))/beta(j);
end
end % tridiag