function [u,it,Vk,Tk] = CGLanczos(A,b,u0,prec,itMax)

u = u0;
r0 = b-A*u;
Anormr0 = sqrt(r0'*A*r0);
Anormr = Anormr0;
r = r0;
p = r;
it = 0;

Vk = [];
Alfa = [];
Beta = [];

while Anormr>prec*Anormr0 && it<itMax
    Vk = [Vk,r/norm(r)];
    alfa = r'*r/(p'*A*p);
    Alfa = [Alfa,alfa];
    u = u + alfa*p;
    r = b-A*u;
    Anormr = sqrt(r'*A*r);
    beta = -r'*A*p/(p'*A*p);
    Beta = [Beta,beta];
    p = r + beta*p;
    it = it+1;
end

diag0 = [1/Alfa(1),1./Alfa(2:end)+Beta(1:end-1)./Alfa(1:end-1)]';
diag1 = ((sqrt(Beta))./Alfa)';
Tk = spdiags([-diag1,diag0,[0;-diag1(1:end-1)]],-1:1,it,it);
