tridiagmatrix <- function ( L, D, U, b )

#*****************************************************************************80
#
## tridiagmatrix solves a tridiagonal linear system.
#
#  Licensing:
#
#    Copyright 2016 James P. Howard, II
#
#    The computer code and data files on this web page are distributed under
#    https://opensource.org/licenses/BSD-2-Clause, the BSD-2-Clause license.
#
#  Modified:
#
#    10 February 2020
#
#  Author:
#
#    Original R code by James Howard;
#    Modifications by John Burkardt.
#
#  Reference:
#
#    James Howard,
#    Computational Methods for Numerical Analysis with R,
#    CRC Press, 2017,
#    ISBN13: 978-1-4987-2363-3.
#
{
  n <- length ( D )
  L <- c ( NA, L )
#
#  Forward sweep.
#
  U[1] <- U[1] / D[1]
  b[1] <- b[1] / D[1]
  for ( i in 2 : ( n - 1 ) )
  {
    U[i] <- U[i] / ( D[i] - L[i] * U[i-1] )
    b[i] <- ( b[i] - L[i] * b[i-1] ) / ( D[i] - L[i] * U[i-1] ) 
  }
  b[n] <- ( b[n] - L[n] * b[n-1] ) / ( D[n] - L[n] * U[n-1] ) 
#
#  Backward sweep.
#
  x <- rep.int ( 0, n )
  x[n] <- b[n]
  for ( i in ( n - 1 ) : 1 )
  {
    x[i] <- b[i] - U[i] * x[i+1]
  }

  return ( x )
}