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