program implicit

implicit none

integer i, j, n, m, INFO
integer, allocatable::IPIV(:)
real ax, bx, dx, dt, t, t1, t2, al
complex f0, st, ima
Complex, allocatable::A(:,:), Ainv(:,:), R(:), Psi(:)
real, allocatable::x(:), y(:), xalt(:), yalt(:), imy(:), rey(:), imyalt(:), reyalt(:)
 
ima=(0.0,1.0)
n=300
m=1000
ax=-20.0
bx=20.0
dx=(bx-ax)/(n-1)
dt=0.02
al = dt/(dx**2)


Allocate(A(n,n), Ainv(n,n), R(n), Psi(n), x(n), y(n), xalt(n), yalt(n),IPIV(n),imy(n), rey(n), imyalt(n), reyalt(n))

do i=1, n
  do j=1,n
    If(i.eq.j)then
    A(i,j)= 1.0+ima*al
    Ainv(i,j)=1.0
    else if((i .eq. (j-1)) .or. (i .eq. (j+1)))then
    A(i,j)=-0.5*ima*al
    Ainv(i,j)=0.0
    else
    A(i,j)=0.0
    Ainv(i,j)=0.0
    end if
  end do
end do

A(1,2)=-ima*al
A(n,n-1)=-ima*al


 
  call CGESV( n, n, A, n, IPIV, Ainv, n, INFO )

  write(*,*) INFO

  call setpag ("DA4P")
  call metafl ("CONS")
  call disini                         
  call THKCRV(5)
 
   call height(50)

do i=1,m

    call messag('t=', 800, 400) 
    
     if(i .gt. 1)then
     call color('BACK')
     call number(t, 2, 900, 400)
     end if

     t=dt*(i-1)
     call color('WHITE')
     call number(t, 2, 900, 400)



      call graf(ax,bx,ax,0.3*bx,-1.2,1.2,-1.2,0.1)
  if(i.eq.1)then
  do j=1,n
  x(j)=ax+dx*(j-1)
  Psi(j)=f0(x(j))
  end do

  else

  do j=1,n 
     R(j) = dt*st(Psi(j)) 
  end do

  Psi = matmul(Ainv,Psi)+matmul(Ainv,R)

  end if

  do j=1,n
    y(j) = abs(Psi(j))**2
    imy(j)=aimag(Psi(j))
    rey(j)=real(Psi(j))
  end do

        Call CPU_TIME(t1)
      
      Do While((t2-t1) .lt. 0.01)
         call CPU_TIME(t2)
      End Do

      If(i .gt. 1) Then         !Löschen der alten Graphen 
      CAll color('BLACK')                       
      Call curve(xalt,yalt,n)
      Call curve(xalt,imyalt,n)
      Call curve(xalt, reyalt,n)
      End If
   
      Call color('WHITE')          
      Call curve(x,y,n)   !Amplitude |Psi|**2 des Solitons
      Call color('GREEN')
      Call curve(x,imy,n) !Imaginärteil
      Call color('RED') 
      Call curve(x,rey,n) !Realteil
      call color('WHITE')
      
      yalt=y
      imyalt=imy
      reyalt=rey
      xalt=x
   
  call endgrf
   
end do

call disfin

end program
function st(x)

implicit none

complex st, x, ima
real g, dt

g=-1.0
ima=(0.0,1.0)

st=-ima*g*abs(x)**2*x

end function
function f0(x)

implicit none

   Real x,x0,a,v
   complex f0, ima


   a=1.0
   x0=0.0
   v=0.0
   ima = (0.0,1.0)
   
   !f0=ima*v+sqrt(a**2-v**2)*tanh(sqrt(a**2-v**2)*(x-x0))
   

   
   f0= a/cosh(a*(x-x0))*exp(-ima*v*x)
   !f0 = a/cosh(a*(x-x0))*exp(-ima*v*x)+a/cosh(a*(x+x0))*exp(ima*v*x)
end function