function Z = adi_single_shift(A, b, tau, k)
n = size(A,1);
v = (A - tau*speye(n)) \ b * sqrt(2*tau);
Z = [v];

% could be improved by computing only once
% a factorization of A - tau*I

for i = 1:k-1
    v = v + (A - tau*speye(n)) \ v * 2*tau;
    Z = [Z v];
end