function x = solve_band5(ssL,sL,d,z)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the linear system B x = z where A is symmetric,
% positive definite and has bandwidth 5.  
% Given the factorization B = LDL', we do two triangular solves:
%    (1)  Solve L u = z for u  (forward substitution)
%    (2)  Solve D v = u
%    (3)  Solve L' x = v for x  (backward substitution)
%
% INPUTS:
%    ssL      Nx1 vector       sub-subdiagonal of L (uses 1st N-2)
%    sL       Nx1 vector       subdiagonal elements of L (uses 1st N-1)
%    d        Nx1 vector       diagonal elements of D
%    z        Nx1 vector       right hand side of equation
%
% OUTPUTS:
%    x        Nx1 vector      solution to B x = z
%
%  Note on indexing of subdiagonal and sub-subdiagonal:
%     index of subdiagonal corresponds to the column e.g., L_21 = sL(1)
%     index of sub-subdiagonal corresponds to column e.g., L_31 = ssL(1)
%
%  Written by Karen Chiswell on 20 March 2006.
%  Follows the exposition in Green & Silverman (1994), sec 2.6.2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = length(z);
%
%  forward solve to get u in L u = z
%
u = z;
u(2) = z(2) - sL(1)*u(1);

for i = 3:N
    u(i) = z(i) - sL(i-1)*u(i-1) - ssL(i-2)*u(i-2);
end   %  forward substitution

%
%  divide by diagonals
%
v = u ./ d;

%
%  backward solve to get x in L' x = v
%
x = v;
x(N-1) = v(N-1) - sL(N-1)*x(N);

for i = N-2:-1:1
    x(i) = v(i) - sL(i)*x(i+1) - ssL(i)*x(i+2);
end  %  backward substitution

% end of solve_band5.m