function x = ConjugateGradient(A, b, x0, tol, maxIters)
    x = x0;
    r = b - A * x;      % initial residual
    p = r;              % search direction set to initial residual
    rsold = r' * r;     % 2-norm squared of the residual

    for i = 1:maxIters
        Ap = A * p;
        alpha = rsold / (p' * Ap);  % correction factor of the solution along direction p
        x = x + alpha * p;  % new approximation
        r = r - alpha * Ap; % new residual (is equal to b-A*x)
        rsnew = r' * r;     % 2-norm squared of the residual

        if sqrt(rsnew) < tol
            break;
        end
        beta = (rsnew / rsold); % this factor guarantees that the new search direction is conjugate to the previous one
        p = r + beta * p;   % new searching direction (is conjugate to all the previous ones)
        rsold = rsnew;
    end
end