function [G, vres] = itfun_U(A, X0, tol, maxiter)
% Fixed point iteration based in matrix U
% On input: X0 = starting approximation
% tol = tolerance for the stopping condition
% maxiter = max number of iterations
% A = 3-dimensional array with block coefficients
% On output: G = approximation of the solution of
% X = A(:,:,1) + A(:,:,1)*X + ... + A(:,:,end)*X^(end-1)
% vres = residuals vector 
k = size(A,3);
m = size(A, 1);
step = 0;
res = tol + 1;
G = X0;
while (res>tol & step<maxiter)
    step = step + 1;
    X1 = A(:,:,k);
    for j = k-1: -1 : 2
        X1 = X1 * G + A(:,:,j);
    end
    R = eye(m) - X1;
    X1 = inv(R)*A(:,:,1);
    res = norm(G-X1, 1);
    if step== 1
        vres = res;
    else
        vres=[vres,res];
    end
    step, res
    G = X1;
end