function [G, vres] = itfun_T(A, X0, tol, maxiter)
% Traditional fixed point iteration
% 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;
R = inv(eye(m)-A(:,:,2));
G = X0;
while (res>tol & step<maxiter)
    step = step + 1;
    X1 = A(:,:,k);
    for j = k-1: -1 : 3
        X1 = X1 * G + A(:,:,j);
    end
    X1 = X1*G*G + A(:,:,1);
    X1 = R*X1;
    res = norm(G-X1, 1);
    if step== 1
        vres = res;
    else
        vres=[vres,res];
    end
    step, res
    G = X1;
end