Short answer: the posted link is right, Maple is wrong.
Long answer: there is a way to determine which of the two proposed formulas cannot be the right answer.
If $\mathcal{L}_X\nabla$ is a tensor, then it must be $\mathcal{C}^{\infty}$-linear in its entries, and in particular in $Z$.
So, let's take $X$, $Y$ and $Z$ three vector fields, and $f$ a smooth function.
We compute:
\begin{align}
\mathcal{L}_X(\nabla_Y(fZ))
&= \mathcal{L}_X((Yf)Z + f\nabla_YZ) \\
&= (X(Yf))Z + (Yf)\mathcal{L}_XZ + (Xf)\nabla_YZ + f\mathcal{L}_X(\nabla_YZ),
\label{eq:1} \tag{1} \\
\\
\nabla_{\mathcal{L}_XY}(fZ)
&= ((\mathcal{L}_XY)f)Z + f \nabla_{\mathcal{L}_XY}Z,
\label{eq:2} \tag{2} \\
\\
\nabla_X(\mathcal{L}_Y(fZ))
&= \nabla_X((Yf)Z + f\mathcal{L}_YZ) \\
&= X(Yf)Z + (Yf)\nabla_XZ + (Xf)\mathcal{L}_YZ + f\nabla_X(\mathcal{L}_YZ),
\label{eq:3} \tag{3} \\
\\
\nabla_Y(\mathcal{L}_X(fZ))
&= \nabla_Y((Xf)Z +f \mathcal{L}_XZ) \\
&= (Y(Xf))Z + (Xf)\nabla_YZ + (Yf)\mathcal{L}_XZ + f\nabla_Y(\mathcal{L}_XZ).
\label{eq:3'} \tag{3'}
\end{align}
Since $(\mathcal{L}_XY)f = X(Yf) - Y(Xf)$, then
$$
\eqref{eq:1}-\eqref{eq:2}-\eqref{eq:3'} = f\mathcal{L}_X(\nabla_YZ) -f\nabla_{\mathcal{L}_XY}Z - f\nabla_Y(\mathcal{L}_XZ),
$$
while
\begin{align}
\eqref{eq:1} - \eqref{eq:2} - \eqref{eq:3}
&= -((\mathcal{L}_XY)f)Z +(Yf)(\mathcal{L}_XZ - \nabla_XZ) + (Xf)(\nabla_YZ-\mathcal{L}_YZ)\\
&\quad + f \mathcal{L}_X(\nabla_YZ) -f\nabla_{\mathcal{L}_XY}Z -f\nabla_X(\mathcal{L}_YZ).
\end{align}
It is clear that $\eqref{eq:1} - \eqref{eq:2} - \eqref{eq:3'}$ is $\mathcal{C}^{\infty}$-linear while $\eqref{eq:1} - \eqref{eq:2} - \eqref{eq:3}$ is not.