function y=phi(tau,A,v)

p=size(v,2)-1;
J=spdiags(ones(p,1),1,p,p);
W=v(:,end:-1:2);
nw=norm(W,inf);
eta = 2^-min(ceil(log2(max(nw,realmin))));
B = [A,eta*W;zeros(p,size(A,2)),J];
ep=flipud(eye(p,1));
v0=[v(:,1); ep/eta];
y=full([speye(size(A)),sparse(length(A),p)]*expm(tau*B)*v0);