	subroutine gradient(psi,psi_gd)
	use constant
	implicit none
	complex(8),intent(in) :: psi(0:nx)
	complex(8),intent(out) :: psi_gd(0:nx)

	complex(8) :: psi1(0:nx),psi2(0:nx),psi3(0:nx),psi4(0:nx)

	if(order.eq.2)  then
		call split_4th(psi,psi1,1)
		call split_4th(psi,psi2,-1)
		psi_gd = (psi1-psi2)/(2.d0*per*eps_para)

	elseif(order.eq.4) then

		call split_4th(psi,psi1,1)
		call split_4th(psi1,psi2,1)
    		psi_gd  = -psi2+8.d0*psi1

		call split_4th(psi,psi1,-1)
		call split_4th(psi1,psi2,-1)
    		psi_gd  =(psi_gd+psi2-8.d0*psi1)/(1.2d1*per*eps_para)

	elseif(order.eq.6) then

		call split_4th(psi,psi1,1)
		call split_4th(psi1,psi2,1)
		call split_4th(psi2,psi3,1)
    		psi_gd  = psi3-9.d0*psi2+4.5d1*psi1

		call split_4th(psi,psi1,-1) ! backward 
		call split_4th(psi1,psi2,-1)
		call split_4th(psi2,psi3,-1)
    		psi_gd  = (psi_gd-4.5d1*psi1+9.d0*psi2-psi3)/(6.d1*per*eps_para)

	elseif(order.eq.8) then

		call split_4th(psi,psi1,1)
		call split_4th(psi1,psi2,1)
		call split_4th(psi2,psi3,1)
		call split_4th(psi3,psi4,1)
		psi_gd = 3.d0*psi4-3.2d1*psi3+1.68d2*psi2-6.72d2*psi1


		call split_4th(psi,psi1,-1) ! backward 
		call split_4th(psi1,psi2,-1)
		call split_4th(psi2,psi3,-1)
		call split_4th(psi3,psi4,-1)

		psi_gd = (psi_gd-(3.d0*psi4-3.2d1*psi3+1.68d2*psi2-6.72d2*psi1))/(-8.4d2*Per*eps_para)

	endif

	end subroutine gradient

	subroutine split(psi_ini,psi_out,flag)
	use constant
	implicit none
	integer(4),intent(in) :: flag
	complex(8),intent(in) :: psi_ini(0:nx)
	complex(8),intent(out) :: psi_out(0:nx)

	integer(4) :: i,j,k
	real(8) :: fr(0:nx),fm(0:nx)
	real(8) :: fr1(0:nx),fm1(0:nx)

	if(flag.eq.1) then
		psi_out = psi_ini
		do k = 1,ntmicro
			do j = 0,nx
				psi_out(j) = exp(-im*eps_para*dtmicro/2.d0*abs(psi_out(j))**2.d0)*psi_out(j)
			enddo

			do j = 0,nx
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,0)
			do j = 0,nx
				psi_out(j) = exp(-im*dtmicro*mu(j))*dcmplx(fr1(j),fm1(j))
			enddo

			do j = 0,nx
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,1)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,1)

			do j = 0,nx
				psi_out(j) = dcmplx(fr1(j),fm1(j))
			enddo

			do j = 0,nx
				psi_out(j) = exp(-im*eps_para*dtmicro/2.d0*abs(psi_out(j))**2.d0)*psi_out(j)
			enddo
		enddo	
	elseif(flag.eq.-1) then

		psi_out = psi_ini
		do k = 1,ntmicro

			do j = 0,nx
				psi_out(j) = exp(im*eps_para*dtmicro/2.d0*abs(psi_out(j))**2.d0)*psi_out(j)
			enddo
			do j = 0,nx
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,0)
			do j = 0,nx
				psi_out(j) = exp(im*dtmicro*mu(j))*dcmplx(fr1(j),fm1(j))
			enddo

			do j = 0,nx
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,1)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,1)

			do j = 0,nx
				psi_out(j) = dcmplx(fr1(j),fm1(j))
			enddo

			do j = 0,nx
				psi_out(j) = exp(im*eps_para*dtmicro/2.d0*abs(psi_out(j))**2.d0)*psi_out(j)
			enddo
		enddo	

	endif

	end subroutine split

	subroutine hermiteForward(psi_phy,psi_spe)
!	
!	Transform from physics to spectral space
!
	use constant
	implicit none
	complex(8),intent(in) :: psi_phy(0:nx)
	complex(8),intent(out) :: psi_spe(0:nx)

	integer(4) :: j
	real(8) :: fr(0:nx),fm(0:nx)
	real(8) :: fr1(0:nx),fm1(0:nx)


	do j = 0,nx
		fr(j) = dreal(psi_phy(j))
		fm(j) = dimag(psi_phy(j))
	enddo
	call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
	call hetrfm1(nx+1,fm,fm1,ax,gxw,0)

	do j = 0,nx
		psi_spe(j) = dcmplx(fr1(j),fm1(j))
	enddo

	end subroutine hermiteForward

	subroutine split_4th(psi_ini,psi_out,flag)
	use constant
	implicit none
	integer(4),intent(in) :: flag
	complex(8),intent(in) :: psi_ini(0:nx)
	complex(8),intent(out) :: psi_out(0:nx)

	integer(4) :: i,j,k
	real(8) :: fr(0:nx),fm(0:nx)
	real(8) :: fr1(0:nx),fm1(0:nx)

	if(flag.eq.1) then
		psi_out = psi_ini
		do k = 1,ntmicro

			do j = 0,nx
				psi_out(j) = exp(-im*eps_para*dtmicro*2.d0*w1*abs(psi_out(j))**2.d0)*psi_out(j)
			enddo

			do j = 0,nx
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,0)
			do j = 0,nx
				psi_out(j) = exp(-im*2.d0*w2*dtmicro*mu(j))*dcmplx(fr1(j),fm1(j))
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,1)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,1)
			do j = 0,nx
				psi_out(j) = dcmplx(fr1(j),fm1(j))
			enddo

			do j = 0,nx
				psi_out(j) = exp(-im*eps_para*dtmicro*2.d0*w3*abs(psi_out(j))**2.d0)*psi_out(j)
			enddo


			do j = 0,nx
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,0)
			do j = 0,nx
				psi_out(j) = exp(-im*2.d0*w4*dtmicro*mu(j))*dcmplx(fr1(j),fm1(j))
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,1)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,1)
			do j = 0,nx
				psi_out(j) = dcmplx(fr1(j),fm1(j))
			enddo

			do j = 0,nx
				psi_out(j) = exp(-im*eps_para*dtmicro*2.d0*w3*abs(psi_out(j))**2.d0)*psi_out(j)
			enddo

			do j = 0,nx
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,0)
			do j = 0,nx
				psi_out(j) = exp(-im*2.d0*w2*dtmicro*mu(j))*dcmplx(fr1(j),fm1(j))
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,1)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,1)
			do j = 0,nx
				psi_out(j) = dcmplx(fr1(j),fm1(j))
			enddo

			do j = 0,nx
				psi_out(j) = exp(-im*eps_para*dtmicro*2.d0*w1*abs(psi_out(j))**2.d0)*psi_out(j)
			enddo

		enddo	


	elseif(flag.eq.-1) then


		psi_out = psi_ini
		do k = 1,ntmicro
			do j = 0,nx
				psi_out(j) = exp(im*eps_para*dtmicro*2.d0*w1*abs(psi_out(j))**2.d0)*psi_out(j)
			enddo

			do j = 0,nx
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,0)
			do j = 0,nx
				psi_out(j) = exp(im*2.d0*w2*dtmicro*mu(j))*dcmplx(fr1(j),fm1(j))
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,1)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,1)
			do j = 0,nx
				psi_out(j) = dcmplx(fr1(j),fm1(j))
			enddo

			do j = 0,nx
				psi_out(j) = exp(im*eps_para*dtmicro*2.d0*w3*abs(psi_out(j))**2.d0)*psi_out(j)
			enddo


			do j = 0,nx
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,0)
			do j = 0,nx
				psi_out(j) = exp(im*2.d0*w4*dtmicro*mu(j))*dcmplx(fr1(j),fm1(j))
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,1)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,1)
			do j = 0,nx
				psi_out(j) = dcmplx(fr1(j),fm1(j))
			enddo

			do j = 0,nx
				psi_out(j) = exp(im*eps_para*dtmicro*2.d0*w3*abs(psi_out(j))**2.d0)*psi_out(j)
			enddo

			do j = 0,nx
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,0)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,0)
			do j = 0,nx
				psi_out(j) = exp(im*2.d0*w2*dtmicro*mu(j))*dcmplx(fr1(j),fm1(j))
				fr(j) = dreal(psi_out(j))
				fm(j) = dimag(psi_out(j))
			enddo
			call hetrfm1(nx+1,fr,fr1,ax,gxw,1)
			call hetrfm1(nx+1,fm,fm1,ax,gxw,1)
			do j = 0,nx
				psi_out(j) = dcmplx(fr1(j),fm1(j))
			enddo

			do j = 0,nx
				psi_out(j) = exp(im*eps_para*dtmicro*2.d0*w1*abs(psi_out(j))**2.d0)*psi_out(j)
			enddo
		enddo	

	endif

	end subroutine split_4th
