How to solve the Schrödinger equation with time dependent potential in 1D or 3D (if it is easier):
$$i\hbar\dfrac{\partial \Psi}{\partial t}(x,t)=\left(-\dfrac{\hbar}{2m}\nabla^2-\frac{e^2}{x+\alpha}-exE(t)\right)\Psi(x,t)$$
where $E(t) = E_0 \exp(-t/\tau^2)sin(\omega_0 t)$ is a Gaussian pulse in time, $\alpha$ is a constant and $e$ is a constant (the electron charge). $\Psi(x,0)$ is hydrogen ground state.
What would it mean to find the solution in a self consistent manner?
On
I'm not sure exactly what problem you are trying to solve (or even if you are still trying to solve it...) but I have used the Crank-Nicolson method to solve something similar. Rather than using a ground-state potential explicitly you might be able to approximate it using a soft Coulomb potential (similar to the second term in the brackets in your question).
Taylor expanding the wavefunction to approximate the time evolution and comparing the result with a time propagator (given by $\hat{T}=\exp\{-(i/\hbar)\hat{H}(t-t_{0})\}$) allows for forward temporal evolution of the wavefunction in the potential to be estimated:
\begin{equation} |\psi(x, t_{0}+t)\rangle =\left[\hat{I}-\frac{i\Delta t}{\hbar}\left(-\frac{\hbar^{2}}{2m_{e}}\frac{\partial^{2}}{\partial x^{2}}+V(x,t)\right)\right]|\psi(x, t_{0})\rangle \end{equation}
$\hat{I}$ is the identity operator, and $\Delta t=(t-t_{0})$. To calculate numerically the second spatial derivative in the Schrödinger equation take the Taylor expansion of both $|\psi(x_{0}+x,t)\rangle$ and $|\psi(x_{0}-x,t)\rangle$ and add them. Rearrange the result for the estimated spatial derivative:
\begin{equation} \frac{\partial^{2}\psi}{\partial x^{2}} = \frac{|\psi(x_{0}+x,t)\rangle - 2|\psi(x,t)\rangle + |\psi(x_{0}-x, t)\rangle}{\Delta x^{2}} \end{equation}
where $\Delta x$ is the numerical spatial step size. You can substitute the spatial derivative in the Schrödinger equation for:
\begin{equation} \psi_{n+1}^{i} = \left[\hat{I}-\frac{i\Delta t}{\hbar}\left(-\frac{\hbar^{2}}{2m_{e}}\frac{\psi_{n}^{i+1} - 2\psi_{n}^{i} + \psi_{n}^{i-1}}{\Delta x^{2}}+V(x,t)\right)\right]\psi_{n}^{i}. \end{equation}
This is the "forward time central step" part of the Crank-Nicolson algorithm, i.e., the forward Euler method (I have simplified the notation so that $|\psi(x_{i},t_{n})\rangle = \psi_{n}^{i}$ here to prevent line wrapping). To get the "backward time central step" term, simply change the sign after $\hat{I}$, which is equivalent to reversing the time propagator (hence "backward time"). The Crank-Nicolson solution to the time-dependent Schrödinger equation is then found by solving:
\begin{equation} \psi_{n+1}^{i} = \left(\hat{I} + \frac{i\Delta t}{2\hbar}\hat{H}\right)^{-1}\left(\hat{I}-\frac{i\Delta t}{2\hbar}\hat{H}\right)\psi_{n}^{i}. \end{equation}
It is a good idea to use sparse matrices if you do intend to solve this numerically as you need a lot of grid points in the time domain and a large spatial window.
This can be quite a cool toy for modelling any time-dependent wavefunction motion for reasonably simple wavefunction evolution, e.g., free-electron wavepacket motion in DC and AC fields, but it can also be used for highly complex light-matter interactions such as high-harmonic generation. For example, the figure below shows a simple single atom model of high-harmonic generation by a very high intensity laser pulse (white line) which is contributing enormously to an electronic potential given by $V(x, t)=V_{\text{soft Coulomb}}(x,t) + V_{\rm{E-Field}}(x,t)$. The colourmap is $\text{log}_{10}(|\psi|^{2})$.
this is not an easy topic to discus, since it would take some writing-up.. you can find in Google the detailed methods. $1D$ is simpler. You can use the split operator time-propagator, where the $x$ and $d/dx$-alike operators are treated differently. When applying the derivative operators, usually one makes use of the Fourier Transform, for which very good software packages are available as open source (FFTW - the Fastest Fourier Transform in the West). Not sure there is any meaning to "consistent manner" in this context . You need to specify the initial wave-function though.