      SUBROUTINE INITAL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NG=400)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/   V(0:NG),   W(0:NG)
     2      /F/  VT(0:NG),  WT(0:NG)
     3      /S/  VX(0:NG),  WX(0:NG), VXX(0:NG), WXX(0:NG),V2W2(0:NG),
     4            X(0:NG)
     5      /C/        XL,        XR,         Q,     ROOT2,      SSET
     6      /I/        IP
C...
C...  LOWER AND UPPER LIMITS OF X
      XL=-30.0D0
      XR= 70.0D0
C...
C...  Q IN CUBIC SCHRODINGER EQUATION
      Q=1.0D0
C...
C...  PRECOMPUTE 2**0.5
      ROOT2=DSQRT(2.0D0)
C...
C...  GRID SPACING
      STEP=(XR-XL)/DFLOAT(NG)
C...
C...  INITIAL REAL, IMAGINARY PARTS (V(X,T), W(X,T)), ABSOLUTE VALUE OF
C...
C...     U(X,0) = U0(X) = (2**0.5)EXP(0.5*I*X)SECH(X), I = (-1)**0.5
C...
C...  OVER XL GE X LE XU
C...
      DO 10 I=0,NG
C...
C...     INITIAL UNIFORM GRID
         X(I)=XL+DFLOAT(I)*STEP
C...
C...     SECH(X)
         SCH=2.0D0/(DEXP(X(I))+DEXP(-X(I)))
C...
C...     REAL PART
         V(I)=ROOT2*DCOS(0.5D0*X(I))*SCH
C...
C...     IMAGINARY PART
         W(I)=ROOT2*DSIN(0.5D0*X(I))*SCH
C...
C...     ABSOLUTE VALUE
         V2W2(I)=DSQRT(V(I)**2+W(I)**2)
10    CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
C...
C...  INITIALIZE COUNTER FOR PLOTTED SOLUTION, CUMULATIVE SUM OF SQUARES
C...  OF ERRORS
      IP=0
      SSET=0.0D0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NG=400)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/   V(0:NG),   W(0:NG)
     2      /F/  VT(0:NG),  WT(0:NG)
     3      /S/  VX(0:NG),  WX(0:NG), VXX(0:NG), WXX(0:NG),V2W2(0:NG),
     4            X(0:NG)
     5      /C/        XL,        XR,         Q,     ROOT2,      SSET
     6      /I/        IP
C...
C...  PDES
C...
C...  V
C...   XX
         NL=1
         NU=1
         CALL DSS044(XL,XR,NG+1,V,VX,VXX,NL,NU)
C...
C...  W
C...   XX
         NL=1
         NU=1
         CALL DSS044(XL,XR,NG+1,W,WX,WXX,NL,NU)
C...
C...  ODES AT THE BOUNDARIES
         VT(0)=0.0D0
         WT(0)=0.0D0
         VT(NG)=0.0D0
         WT(NG)=0.0D0
C...
C...  ODES AT THE INTERIOR POINTS
      DO 10 I=1,NG-1
C...
C...     V**2 + W**2
         V2W2(I)=V(I)**2+W(I)**2
C...
C...     VT
         VT(I)=-WXX(I)-Q*V2W2(I)*W(I)
C...
C...     WT
         WT(I)= VXX(I)+Q*V2W2(I)*V(I)
10    CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NG=400)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/   V(0:NG),   W(0:NG)
     2      /F/  VT(0:NG),  WT(0:NG)
     3      /S/  VX(0:NG),  WX(0:NG), VXX(0:NG), WXX(0:NG),V2W2(0:NG),
     4            X(0:NG)
     5      /C/        XL,        XR,         Q,     ROOT2,      SSET
     6      /I/        IP
C...
C...  INCREMENT THE COUNTER FOR THE PLOTTED OUTPUT
      IP=IP+1
C...
C...  MONITOR PROGRESS OF THE SOLUTION
      WRITE(*,1)NORUN,IP,T
1     FORMAT(' NORUN = ',I2,'  IP = ',I3,'  TIME = ',F8.1)
C...
C...  WRITE THE NUMERICAL AND ANALYTICAL SOLUTIONS AND THE DIFFERENCE
C...  BETWEEN THE TWO IN THE NEIGHBORHOOD OF THE SOLITON PEAK.  COMPUTE
C...  THE SUM OF SQUARES OF ERRORS
      WRITE(NO,2)
2     FORMAT(5X,'I',6X,'T',3X,'X(I)',3X,'X(I)-T',3X,'ABS(UN)'
     +       3X,'ABS(UE)',9X,'DIFF')
      SSE=0.0D0
      DO 3 I=1,NG
         IF(ABS(X(I)-T).LT.2.0D0)THEN
            SRN=DSQRT(V2W2(I))
            SRE=DSQRT(EXACT(X(I),T))
            DIFF=SRN-SRE
            WRITE(NO,4)I,T,X(I),X(I)-T,SRN,SRE,DIFF
4           FORMAT(I6,F7.2,F7.1,F9.3,2F10.5,D13.3)
            SSE=SSE+DIFF**2
         END IF
3     CONTINUE
C...
C...  WRITE THE SUM OF SQUARES OF ERRORS
      WRITE(NO,5)SSE
5     FORMAT(/,' SSE = ',D11.3,//)
C...
C...  PLOT THE SOLUTION
      CALL PLOT1
C...
C...  ACCUMULATE THE SUM OF SQUARES OF ERRORS FOR OUTPUT AT THE END
C...  OF THE RUN
      SSET=SSET+SSE
      IF(IP.EQ.7)THEN
         WRITE(NO,6)SSET
6        FORMAT(//,' CUMULATIVE SSE = ',D11.3,//)
      END IF
      RETURN
      END
      SUBROUTINE PLOT1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NG=400)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/   V(0:NG),   W(0:NG)
     2      /F/  VT(0:NG),  WT(0:NG)
     3      /S/  VX(0:NG),  WX(0:NG), VXX(0:NG), WXX(0:NG),V2W2(0:NG),
     4            X(0:NG)
     5      /C/        XL,        XR,         Q,     ROOT2,      SSET
     6      /I/        IP
C...
C...  WRITE A FILE FOR TOP DRAWER PLOTTING
      IF(IP.EQ.1)THEN
         OPEN(2,FILE='T.TOP',STATUS='NEW')
         WRITE(2,12)
12       FORMAT(' SET LIMITS X FROM -30 TO 70 Y FROM 0. TO 1.5',/,
     1          ' SET FONT DUPLEX')
      END IF
C...
C...  WRITE THE CURRENT SOLUTION
      DO 16 I=1,NG
         SR=SQRT(V2W2(I))
C...
C...     USE THE FOLLOWING STATEMENT TO PLOT THE EXACT SOLUTION
C...     SR=DSQRT(EXACT(X(I),T))
         WRITE(2,11)X(I),SR
11       FORMAT(F10.3,F10.5)
16    CONTINUE
      WRITE(2,13)
13    FORMAT(' JOIN 1')
C...
C...  WRITE THE PLOT LABELS
      IF(IP.EQ.7)THEN
      WRITE(2,14)NORUN
14    FORMAT(
     +' SET FONT DUPLEX'
     +,/,' TITLE 3.0 9.5
     + "Figure ',I2,'a: Solitons of the Cubic Schrodinger Equation"'
     +,/,' TITLE LEFT "         |U| = (V(x,t)**2 + W(x,t)**2)**0.5"'
     +,/,' TITLE BOTTOM "x"'
     +,/,' TITLE 4.5 0.75 "|U| of iU  + U   + Q(|U|**2)U = 0 vs x"'
     +,/,' TITLE 4.5 0.65 "         t     xx"'
     +,/,' TITLE 4.0 0.30 "t = 0, 5,..., 30, Q = 1, NG = 400,
     + Uniform Grid"')
C...
C...  NEXT PLOT
      WRITE(2,15)
15    FORMAT('NEW FRAME')
      END IF
      RETURN
      END
      DOUBLE PRECISION FUNCTION EXACT(X,T)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  THE CUBIC SCHRODINGER EQUATION IS
C...
C...     I*U  + U   + Q*(ABS(U)**2)*U = 0                            (1)
C...        T    XX
C...
C...  IF A COMPLEX FUNCTION IS DEFINED AS
C...
C...     U = V + I*W                                                 (2)
C...
C...  SUBSTITUTION OF EQUATION (2) IN EQUATION (1) GIVES
C...
C...     I*(V  + I*W ) + (V   + I*W  ) + Q*(V**2 + W**2)*(V + I*W) = 0
C...         T      T      XX      XX
C...
C...  OR, SEPARATING INTO REAL AND IMAGINARY PARTS,
C...
C...     V  + W   + Q*(V**2 + W**2)*W = 0                            (3)
C...      T    XX
C...
C...     W  - V   - Q*(V**2 + W**2)*V = 0                            (4)
C...      T    XX
C...
C...  EQUATIONS (3) AND (4) ARE TO BE INTEGRATED SUBJECT TO THE INITIAL
C...  CONDITION
C...
C...     U(X,0) = (2**0.5)*EXP(0.5*I*X)*SECH(X)                      (5)
C...
C...  OR
C...
C...     V(X,0) = (2**0.5)*COS(0.5*X)*SECH(X)                        (6)
C...
C...     W(X,0) = (2**0.5)*SIN(0.5*X)*SECH(X)                        (7)
C...
C...  AND THE BOUNDARY CONDITIONS
C...
C...     U(-INF,T) = U(INF,T) = 0                                 (8)(9)
C...
C...  OR
C...
C...     V(-INF,T) = W(INF,T) = 0                               (10)(11)
C...
C...     V(-INF,T) = W(INF,T) = 0                               (12)(13)
C...
C...  IN THE FOLLOWING CODING, WE TAKE
C...
C...     X = -INF AS X = -30
C...
C...     X =  INF AS X =  70
C...
C...  THE SOLUTION TO THIS PROBLEM IS
C...
C...     U(X,T) = (2**0.5)*EXP(I*(0.5*X + 0.75*T))*SECH(X - T)      (14)
C...
C...  FROM EQUATION (14), WE HAVE
C...
C...     ABS(U(X,T)) = (2**0.5)*SECH(X - T)                         (15)
C...
C...  WHICH IS THE EQUATION FOR THE SOLITON PROGRAMMED BELOW.  NOTE THAT
C...  THE SOLITON TRAVELS FROM LEFT TO RIGHT WITH UNIT VELOCITY AND A
C...  MAXIMUM APLITUDE OF 2**0.5
C...
C...  (ABS(U(X,T)))**2 FROM EQUATION (15)
      EXACT=2.0D0*(2.0D0/(DEXP(X-T)+DEXP(-(X-T))))**2
      RETURN
      END
