function x = tridiagonal_solver ( a, b, c, d )

%*****************************************************************************80
%
%% triadiagonal_solver() solves a tridiagonal linear system.
%
%  Discussion:
%
%    There are N equations in a tridiagonal system.
%
%    Equation 1 is:
%                      b(1) * x(1) + c(1) * x(2)   = d(1)
%    Equation i, for 1 < i < n, is:
%      a(i) * x(i-1) + b(i) * x(i) + c(i) * x(i+1) = d(i)
%    Equation N is:
%      a(n) * x(n-1) + b(n) * x(n)                 = d(n)
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 January 2020
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
%
%  Input:
%
%    real A(N), B(N), C(N), the matrix entries.
%    A(1) and C(N) are not used.
%
%    real D(N,:), the right hand side.
%
%  Output:
%
%    real X(N,:), the solution.
%
  [ m, n ] = size ( b );

  for i = 2 : m
    s = a(i) / b(i-1);
    b(i) = b(i) - s * c(i-1);
    d(i,:) = d(i,:) - s * d(i-1,:);
  end

  x = d;

  for i = m : -1 : 1

    if ( b(i) == 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'tridiagonal_solver - Fatal error!\n' );
      fprintf ( 1, '  b(%d) = 0\n', i );
      error ( 'tridiagonal_solver - Fatal error!' );
    end

    if ( i == m )
      x(i,:) = x(i,:) / b(i);
    else
      x(i,:) = ( x(i,:) - c(i) * x(i+1,:) ) / b(i);
    end

  end

  return
end

