An infinitesimal symmetry of a tangent distribution $D={\rm span}(X_1,\dots,X_n)$ on $n$-dimensional manifold $M$ is a vector field $Y$ such that for every $X\in D$ a Lie bracket $[X,Y]\in D$.
My question is: How do computations look like in a case where $D={\rm span}(X_1:=\partial_z)$, where $M=(x,y,z,t)$? How to find all infinitesimal symmetries of $D$?
You need to be a bit careful: you really mean that for any $X\in\Gamma(D)$ (i.e. $X$ a section of $D$, a linear combination of the vector fields $X_i$ with possibly non-constant coefficients), then $[X,Y] \in\Gamma(D)$. In your case, the general section looks like $$X = f(x,y,z,t)\,\partial_z.$$ Then you should take an arbitrary vector field $$ Y = a(x,y,z,t)\,\partial_x + b(x,y,z,t)\,\partial_y + c(x,y,z,t)\,\partial_z + d(x,y,z,t)\,\partial_t,$$ and find under what conditions $[X,Y]\in\Gamma(D)$, i.e. $[X,Y]\propto \partial_z$. The condition turns out to be that $a, b, d$ don't depend on $z$.