function x = BlockSol(G,H,b)
% G and H are cell arrays produced by BlockTriDiag(A,C)
% b is a N-by-1 cell array that encodes a block vector with m-by-1 entries.
% x is a N-by-1 cell array that encodes a block vector with m-by-1 entries
%    that solves the system Atilde*x = b.
N = length(G);
y = cell(N,1);
x = cell(N,1);

% Forward elimination to solve Gtilde*y = b...
y{1} = G{1}\b{1};
for i=2:N
    y{i} = G{i}\(b{i} - H{i-1}*y{i-1});
end
% Back substitution to solve Gtilde'*x = y...
x{N} = G{N}'\y{N};
for i=N-1:-1:1
    x{i} = G{i}'\(y{i} - H{i}'*x{i+1});
end
