I am not sure exactly what you are asking, but I can give a bit of detail about how we might implement a split operator method. I apologize if this is a orthogonal to exactly what you were asking.
As you well know, time evolution in QM obeys the Schrodinger equation
\begin{equation}
i \hbar \frac{\partial\psi(t)}{\partial t} = \mathcal{H} \psi(t).
\end{equation}
Here, $\mathcal{H}$ is the Hamiltonian of the system in question. Usually, we can write $\mathcal{H}$ as the sum of a kinetic term and a few others, such as an interaction term or a external potential. Here I will focus only on the case of an external potential. In this case we may write
\begin{equation}
\mathcal{H} = \frac{\bar{P}^2}{2m} + V(x, t)
\end{equation}
Time evolution in QM is unitary, so there is some unitary transformation $\hat{U}$ that takes $\psi(x, 0) \rightarrow \psi(x,t)$. We can plug this into the Schrodinger equation to see that:
\begin{align*}
i \hbar \frac{\partial\psi(t)}{\partial t} &= \mathcal{H} \psi(t) \\
i \hbar \frac{\partial\hat{U} \psi(t)}{\partial t} &= \mathcal{H} \hat{U} \psi{(x, 0)}
\end{align*}
Since this is true for an arbitrary $\psi$, it must be true for $\hat{U}$ itself. Setting $\hbar=1$ for ease of notation, we then have the following relation:
\begin{equation}\label{eqn:unitschro}
i \frac{\partial \hat{U}}{\partial t} = \mathcal{H} \hat{U}
\end{equation}
with solution
\begin{equation}
\hat{U} = e^{-i \int_{0}^t \mathcal{H}(t') dt'}
\end{equation}
Thus if we initialize a state to $\psi(x,0)$ we may evolve it by multiplying by the matrix exponential above, or:
\begin{equation}
\psi(x, t) = e^{-i \mathcal{H} t} \psi(x, 0).
\end{equation}
if $\mathcal{H}$ does not depend on $t$.
Now comes actually splitting the operator. Recall that we can write \begin{equation}
\mathcal{H} = \frac{\bar{P}^2}{2m} + V(x, t)
\end{equation}
Recalling that $P \rightarrow -i\nabla$ in the position basis, and that in general $P$ and $V$ will not commute, we may write
\begin{equation}
e^{-i \left(\frac{\bar{P}^2}{2m} + V(x, t)\right) \Delta t} = e^{\frac{i \Delta t}{4m}\nabla^2} e^{-i V \Delta t} e^{\frac{i \Delta t}{4m}\nabla^2} + \mathcal{O}(t^3) \label{eqn:splitoperator}
\end{equation}
where $\Delta t$ is our discretised time step.
The higher order terms on the right come from the fact that the kinetic and potential terms will not commute. This can be seen by expanding the left hand side using the series definition of the exponential, and rewriting the result as a Dyson series. I believe it can also be shown using the Baker-Campell-Hausdorff theorem. In any case, the point is that our original matrix exponential can be written in the above form. This is useful for numerical evolving the wavefunction due to the following observation: In momentum space, the $e^{\frac{i\Delta t}{4m}\nabla^2}$ term is a diagonal matrix which is computationally inexpensive to multiply by, and in position space the same is true for $e^{-iV\Delta t}$. So a rough algorithm to evolve the initial state $\psi(x, 0)$ would be:
Initialize grid
Initialize(psi(x,0))
Initialize (V(x,0))
psi(k) = FFT(\psi(x))
for each timestep:
psi(k) = e^{\frac{i\ Delta t}{4m}k^2} * psi(k)
psi(x) = InverseFFT(psi(k))
psi(x) = e^{-iV\Delta t} * psi(x)
psi(k) = FFT(psi(x))
psi(k) = e^{\frac{i\Delta t}{4m}k^2} * psi(k)
psi(x) = InverseFFT(psi(k))
Essentially, each time step consists of propagating the free particle of the same mass over half a time interval, multiplying by a phase shift that comes from the potential applied over the entire interval, and finishing the step off with a final half time step of free particle propagation. When you have multiple time steps strung together, the final half step of the previous step combines with the first half step of the current one, so you essentially just apply successive full-step propagations alternating between real and momentum space. In general this is not so bad to do since it is limited (usually) by how fact you can do an FFT for whatever grid size you are working with. Very efficient algorithms (FFTW for example) exist for doing this.
