Program Crank_Nichilson
   
   Implicit None
   
   Integer i, j , n, INFO, m	
   Real ax, bx, dx, dt, al,al2, t1, t2, t
   Complex ima, Psi0, f,NOM
   Real, Allocatable::x(:),xalt(:),Psiout(:),Psioutalt(:), phi(:), chi(:), chialt(:), phialt(:), y(:)
   Complex, Allocatable::A(:,:),Ainv(:,:),B(:,:),C(:,:), Psi(:)
   integer,allocatable::IPIV(:)
   
   ima=(0.0,1.0)
   n=200
   m=1000
   ax=-20
   bx=20
   dx=(bx-ax)/(n-1) 
   dt=0.05
   al = dt/(dx**2)
   
   Allocate(A(n,n),Ainv(n,n),B(n,n),C(n,n),Psi(n),x(n),Psiout(n),Psioutalt(n),xalt(n),phi(n),chi(n), phialt(n),chialt(n),IPIV(n),y(n))


   Do i=1,n
      Do j=1,n
         If(i .eq. j) Then
      
         A(i,j) = 0.5*ima*al+(1.0,0.0)
         B(i,j) = -0.5*ima*al+(1.0,0.0)
         Ainv(i,j) = (1.0,0.0)
      Else If( (i .eq. j-1) .Or. (i .eq. j+1))Then
         A(i,j) = -0.25*ima*al
         B(i,j) = 0.25*ima*al
         Ainv(i,j)=(0.0,0.0)
      Else
         A(i,j)=(0.0,0.0)
         B(i,j)=(0.0,0.0)
         Ainv(i,j)=(0.0,0.0)
         end if
      End Do
   End Do
   
   A(1,2)=-0.5*ima*al
   B(1,2)=0.5*ima*al    !Neumann-randbedingung: Ableitung an den Rändern null
   A(n,n-1)=-0.5*ima*al
   B(n,n-1)=0.5*ima*al

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

  write(*,*) INFO


   C= Matmul(Ainv,B)
   
   call setpag("DA4P")
   call metafl("cons")
   call disini
   call height(50)
   call THKCRV(5)
   
   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,10.0,-1.5,1.5,-1.5,0.1)

      If(i .eq. 1) Then
      Do j=1,n
         x(j)   = ax+dx*(j-1)
         
         Psi(j) = Psi0(x(j))         !Anfangsverteilung
         !y(j)   = 0.005*x(j)**2    
         !y(j)   = 0.5*sin(0.3*x(j))**2 
          !y(j)  = -0.008*x(j)**2+0.00003*x(j)**4+0.3 
         phi(j) = real(Psi(j))
         chi(j) = aimag(Psi(j))
         Psiout(j)=abs(Psi(j))**2              
      end do
      
      
   Else
      
      Do j=1,n
         Psi(j) = f(Psi(j),x(j),dt)*Psi(j)    !1/3 schritte
      End Do


      Psi=matmul(C,Psi)         !2/3 schritte

      do j=1,n
	   Psi(j) = f(Psi(j),x(j),dt)*Psi(j)
           phi(j)=real(Psi(j))             !letzte schritt
           chi(j)=aimag(Psi(j))
           Psiout(j)=abs(Psi(j))**2  
      end do
     
      
      
      
      end if
      
      
      
      
      
      
      Call CPU_TIME(t1)
      
      Do While((t2-t1) .lt. 0.01)
         call CPU_TIME(t2)
      End Do
      
      If(i .gt. 1) Then
      CALL color('BLACK')
      call curve(xalt,Psioutalt,n) 
      call curve(xalt,phialt,n)
      call curve(xalt,chialt,n)
      End If
      

      CALL color('WHITE')
      call curve(x,Psiout,n)
      call color('RED')
      call curve(x,phi,n)
      call color('GREEN')
      call curve(x,chi,n)
      call color('WHITE')
      call curve(x,y,n)
      
      xalt=x
      Psioutalt=Psiout
      phialt= phi
      chialt= chi
      
      call endgrf
      
   End Do
   
   CALL disfin
   
   
   
   
   
End Program
Function f(psi,x,dt)
   Implicit None
   Real dt, g, x
   Complex f, ima,psi
   
   g=-1.0
   ima=(0.0,1.0)
   

   f=EXP(-ima*dt*g*0.5*Abs(psi)**2)   
   !f=exp(-ima*dt*0.5*(-0.008*x**2+0.00003*x**4+0.5+g*abs(psi)**2))
   !f=EXP(-ima*dt*0.5*(0.005*x**2+g*Abs(psi)**2))
End Function
Function Psi0(x)
   Implicit None
   Complex Psi0, ima
   Real x,x0,v,al, a
   
   v=1.0
   ima=(0.0,1.0)
   a=1.0
   x0=0.0

   !if(x.lt.0.0)then
   !Psi0=ima*v+ tanh(x+x0)
   !else if (x.gt. 0.0)then
   !Psi0=ima*v-tanh(x-x0)
   !end if
   
   Psi0 =a/Cosh(a*(x-x0))*Exp(-ima*v*x)   !bright soliton
   !Psi0 =a/Cosh(a*(x-x0))*Exp(-ima*v*x)+a/Cosh(a*(x+x0))*Exp(ima*v*x)
   !Psi0=ima*v+tanh(x-x0)
  
End Function