Program explicit_soliton_ds
   
   Implicit None
   

   Integer i, j, n, m
   Real ax, bx, dx, dt, al,g, be, t1, t2, t
   Real, Allocatable::x(:),phi(:),chi(:),xalt(:),phialt(:),chialt(:),u(:),ualt(:)
   complex ima,Psi0
   complex, allocatable::Psi(:,:)

   n= 200
   m=1000
   ax=-18.0     
   bx=18.0
   dx= (bx-ax)/(n-1)
   dt=0.01
   ima=(0.0,1.0)
   g=1.0
   al=dt/(dx**2)
   be=g*dt
   
   
   
   Allocate(x(n),phi(n),chi(n),Psi(m,n),u(n),ualt(n),xalt(n),phialt(n),chialt(n))
   
   Do i=1,n
      x(i)= ax+dx*(i-1) 
      Psi(1,i)=Psi0(x(i))     !Anfangasbedingung
   End Do
   

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

 
   call height(50)
   
   Do i=2,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.5*bx,-1.2,1.2,-1.2,0.1)
      
      !*******************************************************************************************************

    do j=2,n-1
      if(i .eq. 2)then 
              Psi(i,j)=Psi(1,j)+0.5*(ima*al*Psi(i-1,j+1)-2.0*ima*al*Psi(i-1,j)+1.0*ima*al*Psi(i-1,j-1)-2.0*ima*be*abs(Psi(i-1,j))**2*Psi(i-1,j))
      else 
      Psi(i,j)=Psi(i-2,j)+1.0*ima*al*Psi(i-1,j+1)-2.0*ima*al*Psi(i-1,j)+1.0*ima*al*Psi(i-1,j-1)-2.0*ima*be*abs(Psi(i-1,j))**2*Psi(i-1,j)
      end if
    end do

    
      Psi(i,n)=Psi(i,n-2)         !Neumann-Randbedingung: Ableitung an den Rändern null
      Psi(i,1)=Psi(i,3)
      
      
      Do j=1,n
        chi(j)=aimag(Psi(i,j))       !Ausgabe
        phi(j)=real(Psi(i,j))
        u(j) = abs(Psi(i,j))**2 
      End Do
      
      !*********************************************************************************************************
      
      
      !****************************Warte-Schleife für die graphische Darstellung*******************************************
      Call CPU_TIME(t1)
      
      Do While((t2-t1) .lt. 0.01)
         call CPU_TIME(t2)
      End Do
      !********************************************************************************************************************
      
      !***************************Graphen mit Dislin*********************************************************************
            If(i .gt. 1) Then         !Löschen der alten Graphen 
      CAll color('BLACK')                       
      Call curve(xalt,ualt,n)
      Call curve(xalt,chialt,n)
      Call curve(xalt, phialt,n)
      End If
   
      Call color('WHITE')          
      Call curve(x,u,n)   !Amplitude |Psi|**2 des Solitons
      Call color('GREEN')
      Call curve(x,chi,n) !Imaginärteil
      Call color('RED') 
      Call curve(x,phi,n) !Realteil
      call color('WHITE')
      
      ualt=u
      phialt=phi
      chialt=chi
      xalt=x
      
      !***************************************************************************************************************
     call endgrf()

      
  End Do

CALL disfin() 
   
End Program
Function Psi0(x)
   Implicit None
   Real x,x0,a,v
   complex Psi0, ima

   a=1.0
   x0=0.0
   v=0.0
   ima = (0.0,1.0)
   
   Psi0=ima*v+sqrt(a**2-v**2)*tanh(sqrt(a**2-v**2)*(x-x0))
   
   !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)
   !else
   !Psi0=0.0
   !end if

   
   !Psi0= a/cosh(a*(x-x0))*exp(-ima*v*x)
   !Psi0 = a/cosh(a*(x-x0))*exp(-ima*v*x)+a/cosh(a*(x+x0))*exp(ima*v*x)
End Function