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

%*****************************************************************************80
%
%% tridisolve() solves a tridiagonal system of equations.
%
%  Discussion:
%
%    x = TRIDISOLVE(a,b,c,d) solves the system of linear equations
%    b(1)*x(1) + c(1)*x(2) = d(1),
%    a(j-1)*x(j-1) + b(j)*x(j) + c(j)*x(j+1) = d(j), j = 2:n-1,
%    a(n-1)*x(n-1) + b(n)*x(n) = d(n).
%
%    The algorithm does not use pivoting, so the results might
%    be inaccurate if abs(b) is much smaller than abs(a)+abs(c).
%    More robust, but slower, alternatives with pivoting are:
%      x = T\d where T = diag(a,-1) + diag(b,0) + diag(c,1)
%      x = S\d where S = spdiags([[a; 0] b [0; c]],[-1 0 1],n,n)
%
%  Licensing:
%
%    Copyright 2014 The MathWorks, Inc.
%
%  Author:
%
%    Cleve Moler
%
%  Reference:
%
%    Cleve Moler,
%    Numerical Computing with MATLAB,
%    SIAM, 2004,
%    ISBN13: 978-0-898716-60-3,
%    LC: QA297.M625,
%    ebook: https://www.mathworks.com/moler/chapters.html
%
x = d;
n = length(x);

for j = 1:n-1
   mu = a(j)/b(j);
   b(j+1) = b(j+1) - mu*c(j);
   x(j+1) = x(j+1) - mu*x(j);
end

x(n) = x(n)/b(n);
for j = n-1:-1:1
   x(j) = (x(j)-c(j)*x(j+1))/b(j);
end
