for k = 1:kmax
q0 = copy(q1)
q1 = r/beta
r = A * q1
alpha = dot(q1, r)
T[k,k] = alpha
if k > 1
T[k-1,k] = beta
T[k,k-1] = beta
end
r = r - alpha*q1 - beta*q0
beta = norm(r)
end
This is much more computationally efficient than Arnoldi.