First I will give the relevant information for my question, and then I'll ask the question.
$\large{\textrm{Background}}$
For evolving the nonlinear Schrodinger equation (NLS), one typically uses [a variant of] the Split-Operator, in which one split the nonlinear and linear components of the differential equation and treats them seperately, one after another. The version of the NLS I'm concerned with is
$$\frac{\partial \psi}{\partial t}=\frac{i}{2}\frac{\partial^2 \psi}{\partial x^2}+i|\psi|^2\psi=[\hat{D}+\hat{N}(x,t)]\psi$$
where we have defined
$$\psi=\psi(x,t),\,\,\hat{D}=\frac{i}{2}\frac{\partial^2}{\partial x^2},\,\,\hat{N}(x,t)=i|\psi|^2$$
As outlined quite succinctly in the previously cited wiki page, the general evolution $\psi(x,t)\to\psi(x,t+\Delta t)$ for sufficiently small step size $\Delta t$ is given by
$$\psi(x,t+\Delta t)=e^{\Delta t(\hat{D}+\hat{N})}\psi(x,t)$$
To actually evolve this with a computer, we can split/factor the exponential operator in the line above to a obtain an algorithm that is correct to order $(\Delta t)^2$.
$$\psi(x,t+\Delta t)\approx e^{\Delta t(\hat{D})}e^{\Delta t(\hat{N})}\psi(x,t)$$
Realizing that the second derivative in the $\hat{D}$ operator turns into a multiplication by $-k^2$ in Fourier-space, we can perform this evolution on a computer by calculating, in a discretized fashion of course,
$$\psi(x,t+\Delta t)\approx\mathcal{F}^{-1}\left(e^{-\frac{i}{2}\Delta t (k^2)}\mathcal{F}\left(e^{i\Delta t |\psi^2|}\psi(x,t)\right)[k,t]\right)[x,t]$$
This can be interpreted as follows: we first "evolve" the nonlinear part of the equation, and then we take a time step $\Delta t$ by the "evolution" of the linear part of the equation (in Fourier space). With this scheme, one can start with an initial wave-form/function $\psi(x,t_0)$ and evolve forward in time.
$\large{\textrm{Question}}$
How does one use higher order algorithms on the NLS? By using the Baker-Campell-Hausdorff formula for factorizing exponentials of operators, one can come up with a simple third order algorithm,
$$\psi(x,t+\Delta t)\approx e^{\frac{\Delta t}{2}\hat{D}}e^{\Delta t(\hat{N})}e^{\frac{\Delta t}{2}\hat{D}}\psi(x,t)$$
However, if we try to attack this in the same way we did before, we now have to evaluate $\hat{N}(x,t+\frac{\Delta t}{2})$, which involves evaluating $\psi(x,t+\frac{\Delta t}{2})$, which is the function we are trying to evolve! More precisely, if we only start off with an initial wave-form/function $\psi(x,t)$, then we don't have access to $\psi(x,t+\frac{\Delta t}{2})$. How can I resolve this problem.
Edit: After writing this I realized that by evaluating $e^{\frac{\Delta t}{2}\hat{D}}\psi(x,t)$, it's possible we are obtaining an approximation for $\psi(x,t+\frac{\Delta t}{2})$, which would answer my question affirmatively. However I'm going to keep this question up in case I'm incorrect.
