function [x] = conjgrad(A, b)
    [n,~] = size(A);
    x = zeros(n,1);
    r = b - A * x;
    p = r;
    rsold = r' * r;

    for i = 1:length(b)
        Ap = A * p;
        alpha = rsold / (p' * Ap);
        x = x + alpha * p;
        r = r - alpha * Ap;
        rsnew = r' * r;
        if sqrt(rsnew) < 1e-8
              break;
        end
        p = r + (rsnew / rsold) * p;
        rsold = rsnew;
    end
end