First, determine the HH transformation $Q_1$ such that $Q_1^Tz=\beta e_1$ (where $\beta=\pm\|z\|_2$ depending on the implementation of the HH transformation; you can multiply $Q_1$ by $-1$ if you are not happy with the sign of $\beta$) and set $B:=Q_1^TAQ_1$. Then apply the usual "Hessenbergization" to $B$ instead of $A$: $Q_{n-1}^T\cdots Q_2^TBQ_2\cdots Q_{n-1}=H$. The matrix $Q:=Q_1\cdots Q_{n-1}$ then satisfies $Q^Tz=\beta e_1$ and $Q^TAQ=H$.
Here is a (very simplified but working) MATLAB code:
function [Q, H] = hess2(A, z)
assert(ismatrix(A) && size(A,1) == size(A,2),...
'A must be a square matrix.');
n = size(A, 1);
assert(all(size(z) == [n, 1]),...
'Z must be a %d x 1 vector.', n);
Q = eye(n);
H = A;
for k = 1 : n - 1
if k == 1
v = hhmat(z);
else
v = hhmat(H(k:end,k-1));
end
V = eye(n - k + 1) - 2 * v * v';
W = blkdiag(eye(k - 1), V);
H = W * H * W;
Q = Q * W;
end
function v = hhmat(x)
% Compute the vector V so that H := I - 2 * V * V' is a Householder
% transformation satisfying H * V = s * e_1, where s is + or - 1.
v = x;
alpha = -norm(v);
if (v(1) < 0) alpha = -alpha; end
v(1) = v(1) - alpha;
v = v / norm(v);
Alternatively (as a better option), you can first compute $B:=Q_1^TAQ_1$ and then use the MATLAB's hess function to compute the Hessenberg form of $B$ (and compute the product of the two orthogonal matrices $Q_1$ and that obtained from hess).
The reason why this works is that the HH matrices involved in the "Hessenbegization" have the form
$$
\pmatrix{I_k&0\\0&\tilde{Q}},
$$
where $I_k$ is a $k\times k$ identity and $k\geq 1$. So in
$$
Q^Tz=Q_{n-1}^T\cdots Q_2^TQ_1^T=Q_{n-1}^T\cdots Q_2^T(\beta e_1),
$$
they act only on the zero components of the vector $Q_1^Tz=\beta e_1$.