function [x,relerr] = iSkewSolve(b,c)
% Solves Ax = b where b is n-by-1 and
% A = eye(n,n) + diag(-c(1:n-1),-1) + diag(c(1:n-1),1)
n = length(b);
% Compute the factorization A = LU. L will be unit lower bidiagonal and
% U will be upper bidiagonal.
% d(1:n) will house the diagonal of U.
% c(1:n-1) is U's superdiagonal
% m(1:n-1) will house L's subdiagonal
% Derive algorithm by comparing entries in
%
% 1 c(1) 0 0 1 0 0 0 d(1) c(1) 0 0
% -c(1) 1 c(2) 0 m(1) 1 0 0 0 d(2) 0 0
% 0 -c(2) 1 c(3) = 0 m(2) 1 0 0 0 d(3) c(3)
% 0 0 -c(3) 1 0 0 m(3) 1 0 0 0 d(4)

d(1) = 1;
for k=1:n-1
m(k) = -c(k)/d(k);
d(k+1) = 1-m(k)*c(k);
end
% Now overwrite b with the solution to Ly = b...
for k=2:n
b(k) = b(k) - m(k-1)*b(k-1);
end
% Next, overwrite b with the solution to Ux = b.
b(n) = b(n)/d(n);
for k=n-1:-1:1
b(k) = (b(k)-c(k)*b(k+1))/d(k);
end
x = b;

% Error estimate...
% We solve (A+E)xhat = b.
% relerr in xhat is appx (norm(E)/norm(A))*cond(A) = norm(E)*norm(inv(A))
% where A = I + S.
% Thru experimentation we conclude that norm inv(A) is
% approximately 1
normAinv = 1;
% Using GVL (4.2.1) and (4.2.3) we use
normE = eps*(1 + max(abs(c))^2);
relerr = normE*normAinv;