To simplify this problem to a large extent, we constructed a wave function using the following equation.
\begin{equation}
\psi = A\exp(\frac{-(x-x_o)^2}{\sigma^2})\exp(ikx)
\end{equation}
$A$ is a constant such that $\int \psi^*\psi dx = 1$.
\begin{equation}
i\hbar\frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m}\frac{\partial^2 \psi}{\partial x^2}+V\psi
\end{equation}
In reduced form,
\begin{equation}
\frac{\partial \psi}{\partial t} = \frac{i}{2}\frac{\partial^2 \psi}{\partial x^2}-iV\psi
\end{equation}
By Taylor's expansion, we can approximate $\frac{\partial \psi}{\partial t}$ and $\frac{\partial^2 \psi}{\partial x^2}$.
\begin{equation}\label{eq:1a}
\frac{\partial \psi}{\partial t} = \frac{\psi(x,t+\Delta t)-\psi(x,t)}{\Delta t}
\end{equation}
\begin{equation}\label{eq:1}
\psi(x+\Delta x,t) \approx \psi(x,t)+\frac{\partial \psi}{\partial x}\Delta x + \frac{1}{2!}\frac{\partial^2 \psi}{\partial x^2}\Delta^2 x
\end{equation}
\begin{equation}\label{eq:2}
\psi(x-\Delta x,t) \approx \psi(x,t)-\frac{\partial \psi}{\partial x}\Delta x + \frac{1}{2!}\frac{\partial^2 \psi}{\partial x^2}\Delta^2 x
\end{equation}
Therefore, we will obtain,
\begin{equation}\label{eq:2a}
\frac{\partial^2 \psi}{\partial x^2} = \frac{\psi(x+\Delta x,t)-2\psi(x,t)+\psi(x-\Delta x,t)}{\Delta^2 x}
\end{equation}
A wave function consists of real and imaginary parts.
\begin{equation}
\frac{\partial (\psi_r+i\psi_i)}{\partial t} = \frac{i}{2}\frac{\partial^2 (\psi_r+i\psi_i)}{\partial x^2}-iV(\psi_r+i\psi_i)
\end{equation}
So, we have obtained,
\begin{equation}\label{eq:3}
\frac{\partial \psi_r}{\partial t} = \frac{-1}{2}\frac{\partial^2 \psi_i}{\partial x^2}+V\psi_i
\end{equation}
\begin{equation}\label{eq:4}
\frac{\partial \psi_i}{\partial t} = \frac{1}{2}\frac{\partial^2 \psi_r}{\partial x^2}-V\psi_r
\end{equation}
With Equations \ref{eq:1a} and \ref{eq:2a}, we can rewrite Equations \ref{eq:3} and \ref{eq:4},
\begin{equation}\label{eq:5}
\psi_r(x,t+\Delta t) = \psi_r(x,t)-\frac{\Delta t}{2\Delta^2 x}(\psi_i(x+\Delta x,t)-2\psi_i(x,t)+\psi_i(x-\Delta x,t))+V\psi_i(x,t)\Delta t
\end{equation}
\begin{equation}\label{eq:6}
\psi_i(x,t+\Delta t) = \psi_i(x,t)+\frac{\Delta t}{2\Delta^2 x}(\psi_r(x+\Delta x,t)-2\psi_r(x,t)+\psi_r(x-\Delta x,t))-V\psi_r(x,t)\Delta t
\end{equation}
Also it is important to normalize $\psi_r$ and $\psi_i$ such that $\int \psi^*\psi dx = 1$. With Equations \ref{eq:5} and \ref{eq:6}, we are ready to write a program to solve this problem. To ensure stability, $\Delta x$ has to be always greater than $\Delta t$.