function Z = adi_single_shift(A, b, tau, k)
n = size(A, 1);

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