Here is an observation which does not answer the question, but does at least tell where not to look for examples. It is shown that for compact Riemannian manifolds with non-positive Ricci curvature there are no interesting solutions to the stated equations, so to find interesting solutions one must look: a. in positive curvature, b. at indefinite metrics, c. on noncompact manifolds, or d. in infinite dimensions.
Suppose $M$ is compact and $g$ is Riemannian. It is claimed that if $g$ has non-positive Ricci curvature, then either the Ricci curvature of $g$ is somewhere negative and both $X$ and $Y$ are $0$, or $g$ is Ricci flat and both $X$ and $Y$ are parallel. To begin with, do not suppose anything about the Ricci curvature of $g$. Let $D$ be its Levi-Civita connection and raise and lower indices using $g_{ij}$ and the inverse symmetric bivector $g^{ij}$. In what follows I use the abstract index notations, so square brackets (resp. parentheses) denote anti-symmetrization (resp. symmetrization) over the enclosed indices.
For any metric $(L_{X}g)_{ij} = 2D_{(i}X_{j)}$ and
\begin{equation}\label{e1}
(L_{X}^{2}g)_{ij} = X^{p}D_{p}(L_{X}g)_{ij} + (D_{i}X^{p})(L_{X}g)_{pj} + (D_{j}X^{p})(L_{X}g)_{ip}.
\end{equation}
Trace this to obtain
\begin{equation}\label{e2}
g^{ij}(L_{X}^{2}g)_{ij} = 2X^{p}D_{p}D^{q}X_{q} + 4D^{(p}X^{q)}D_{(p}X_{q)}.
\end{equation}
By assumption this equals $g^{ij}(L_{Y}g)_{ij} = 2D^{p}Y_{p}$. Since by assumption $M$ is compact and without boundary, integration by parts yields
\begin{align}\label{e3}
0 = 4\int_{M}D^{(p}X^{q)}D_{(p}X_{q)} - 2\int_{M}(D_{p}X^{p})^{2}.
\end{align}
In general, this yields no obvious conclusions. Go back to the equation preceeding the integration and commute derivatives to obtain
\begin{align}\label{e4}
g^{ij}(L_{X}^{2}g)_{ij} = 2X^{p}D^{q}D_{p}X_{q} + 4D^{(p}X^{q)}D_{(p}X_{q)} -2R_{pq}X^{p}X^{q},
\end{align}
in which $R_{ij}$ is the Ricci curvature of $g_{ij}$. Now integrating by parts yields
\begin{align}\label{e5}
\begin{split}0 &= \int_{M}\left(4D^{(p}X^{q)}D_{(p}X_{q)} - 2D^{q}X^{p}D_{p}X_{q} - 2R_{pq}X^{p}X^{q}\right)\\
& = 2\int_{M}\left(D^{(p}X^{q)}D_{(p}X_{q)} + D^{[p}X^{q]}D_{[p}X_{q]} -R_{pq}X^{p}X^{q}\right).
\end{split}
\end{align}
Because $g$ is Riemannian, if the Ricci curvature is non-positive this implies that $D_{(i}X_{j)} = 0$ and $D_{[i}X_{j]} = 0$ so that $D_{i}X_{j} = 0$ and $X$ is parallel. In the original equation this implies $Y$ is Killing. By the original Bochner argument, if the Ricci curvature is non-positive and somewhere negative then there is no non-zero Killing field, so the only possibility for non-trivial solutions is that $g$ be Ricci flat, in which case the Bochner argument forces $Y$ to be parallel.
