function [V, H, beta] = arnoldi(A, b, n)
m = length(b);
V = zeros(m, n+1);
H = zeros(n+1, n);

beta = norm(b);
V(:, 1) = b / beta;
for j = 1:n
    w = A * V(:,j);  % continuation vector
    for i = 1:j
        H(i,j) = V(:,i)' * w;
        w = w - V(:,i) * H(i,j);
    end
    H(j+1,j) = norm(w);
    V(:,j+1) = w / H(j+1,j);
end