# Numerical solution to Time-dependent Schrodinger equation with time dependent hamiltonian

Currently I am facing the problem to solve numerically the following equation for a double well harmonic potential:

$iℏ\frac{\partial}{\partial t}ψ(x,t)=−\frac{ℏ}{2m}\frac{\partial ^2}{∂x^2}ψ(x,t)+V(x,t)ψ(x,t)$

where

$V(x,t)=\frac{mω^2_0}{2}[a(t)x^4−b(t)x^2]$

Where both $a(t)$ and $b(t)$ are of the form $\frac{1}{2*\sqrt{const_1}}(const_1-const_2 \cos(\omega t))$. I started to use the split-step Fourier method algorithm and here is my code:

import numpy as np
import matplotlib.pyplot as plt

tfinal=10000#total simulation time
dt=0.001#time step
M=np.round(tfinal/dt)
time_steps=M.astype(int)# number of time steps

b=10.# final x
a=-10.# initial x
N=512# num of pieces

dx=(b-a)/N # x-space step size
dk=2*np.pi/(b-a)
n=np.linspace(-N/2,N/2 -1,N) # this is gonna be my n's in x_n=a+dx*n

x=n*dx
k=n*dk
t=0

alpha=0.0005
beta=0.0001
acoef=(alpha-beta*np.cos(t))/(2*alpha**0.5)
bcoef=(alpha-beta*np.cos(t))**0.5 / (2*alpha**0.5)
V=0.5*(acoef*x**4 -bcoef*x**2)

u0=np.exp(-(x+4.72)**2 /2)*(1/np.pi)**0.25
u=u0

for i in xrange(1,time_steps):
u=np.exp(-1j*dt*V/2)*u
c=np.fft.fftshift(np.fft.fft(u))
c=np.exp(-1j*dt*k**2 /2)*c
u=np.fft.ifft(np.fft.fftshift(c))
u=np.exp(-1j*dt*V/2)*u
t=t+dt

plt.plot(x,u*np.conjugate(u),'r--')
plt.xlim(-10,10)
plt.ylim(0,0.6)
plt.show
`

However the problem that arises it that it takes a very long time to execute in python. I am not sure what am I doing wrong since I am new in python or if I am implementing correctly the algorithm. I would appreciate any insight or suggestion.