%   Thanks   Artem Kultyshev
%[A,x-solutions,residuals,number of iterations] = cg(matrix dim,max number
%of iterations)
function [A,x1,r1,iter] = cg(n,maxit)  
    A = [];
    b = []';
    x0 = [];
    %building initial matrices
    for i = 1:n
        x0(i) = 0;
        b(i) = 1;
        for j = 1:n
            A(i,j) = 1/(i+j-1);
        end
    end
    
    b = b';
    x0 = x0';
    r0 = A*x0 - b;
    p0 = -r0;
    iter = 0;
    for i = 1:maxit
        iter = iter + 1;
        a0 = r0'*r0/(p0'*A*p0);
        x1 = x0 + a0*p0;
        r1 = r0 + a0*A*p0;
        b1 = r1'*r1/(r0'*r0);
        p1 = -r1+b1*p0;
        if (sum(abs(r1)>10^-6) == 0)
            break;  
        end 
        x0 = x1;
        r0 = r1;
        p0 = p1;
    end
end 

    
