function [Lsub,Udiag,Usuper] = trilu_decomp(Asub,Adiag,Asuper)

%  [Lsub,Udiag,Usuper] = trilu_decomp(Asub,Adiag,Asuper)
%
%  Compute LU decomposition of tridiagonal matrix A, without pivoting.
%
%  Input variables:
%    Asub contains subdiagonal of A. Asub(i) = A(i+1,i)
%    Adiag contains main diagonal of A. Adiag(i) = A(i,i)
%    Asuper contains superdiagonal of A. Asuper(i) = A(i,i+1)
%
%  Output variables:
%    Lsub contains subdiagonal of L. L(i) = L(i+1,i).
%    Udiag contains main diagonal of U. Udiag(i) = U(i,i)
%    Usuper contains superdiagonal of U. Usuper(i) = U(i,i+1)

  n = length(Adiag);
  for i=1:n-1
    ip1 = i+1;
    m_i = Asub(i) / Adiag(i);  %  Multiplier.
    Adiag(ip1) = Adiag(ip1) - m_i*Asuper(i);
    Lsub(i) = m_i;  %  Store multiplier in lower subdiagonal of L.
  end

  %  Store modified upper triangular components of A.

  Udiag = Adiag;
  Usuper = Asuper;