C-----------------------------------------------------------------------
c     This is a modification of a problem from Ismail and Taha (2001)
c     It is a nonlinear Schrodinger problem involving two
c     complex solution components that we separate into four
c     real functions as follows:

c     u1 --- real part of the first PDE component
c     u2 --- imaginary part of the first PDE component
c     u3 --- real part of the second PDE component
c     u4 --- imaginary part of the second PDE component 

c     See below for the defintions for the PDE, boundary conditions,
c     and initial conditions. The spatial domain is [-30, 90] and
c     t_0 = 0. It does not have a closed form solution.
c
c-----------------------------------------------------------------------
      subroutine derivf(t, x, u, ux, uxx, dfdu, dfdux, dfduxx, npde)
c-----------------------------------------------------------------------
c Purpose:
c       This subroutine is used to define the information about the 
c       PDE required to form the analytic Jacobian matrix for the DAE
c       or ODE system. Assuming the PDE is of the form
c                        ut = f(t, x, u, ux, uxx)
c       this routine returns the Jacobians d(f)/d(u), d(f)/d(ux), and
c       d(f)/d(uxx).
c
c-----------------------------------------------------------------------
c Subroutine Parameters:
c Input:
        integer                 npde
c                               The number of PDEs in the system.
c
        double precision        t
c                               The current time coordinate.
c
        double precision        x
c                               The current spatial coordinate.
c
        double precision        u(npde)
c                               u(1:npde) is the approximation of the
c                               solution at the point (t,x).
c
        double precision        ux(npde)
c                               ux(1:npde) is the approximation of the
c                               spatial derivative of the solution at
c                               the point (t,x).
c
        double precision        uxx(npde)
c                               uxx(1:npde) is the approximation of the
c                               second spatial derivative of the 
c                               solution at the point (t,x).
c
c Output:
        double precision        dfdu(npde,npde)
c                               dfdu(i,j) is the partial derivative
c                               of the i-th component of the vector f
c                               with respect to the j-th component
c                               of the unknown function u.
c
        double precision        dfdux(npde,npde)
c                               dfdux(i,j) is the partial derivative
c                               of the i-th component of the vector f
c                               with respect to the j-th component
c                               of the spatial derivative of the 
c                               unknown function u.
c
        double precision        dfduxx(npde,npde)
c                               dfduxx(i,j) is the partial derivative
c                               of the i-th component of the vector f
c                               with respect to the j-th component
c                               of the second spatial derivative of the 
c                               unknown function u.
c-----------------------------------------------------------------------

        double precision        epsi, epsi2
c-----------------------------------------------------------------------
c
c     Assign dfdu(1:npde,1:npde), dfdux(1:npde,1:npde), and
c     dfduxx(1:npde,1:npde) according to the right hand side of the PDE 
c     in terms of u(1:npde), ux(1:npde), uxx(1:npde).
c
      epsi = 2.d0/3.d0
      epsi2 = 2.d0 * epsi

      dfdu(1,1) = - 2.d0 * u(1) * u(2)
      dfdu(1,2) = - u(1) * u(1) - 3.d0 * u(2) * u(2) - epsi 
     &            * (u(3) * u(3) + u(4) * u(4))
      dfdu(1,3) = - epsi2 * u(3) * u(2)
      dfdu(1,4) = - epsi2 * u(4) * u(2)
      dfdu(2,1) = 3.d0 * u(1) * u(1) + u(2) * u(2) + epsi
     &            * (u(3) * u(3) + u(4) * u(4))
      dfdu(2,2) = 2.d0 * u(2) * u(1)
      dfdu(2,3) = epsi2 * u(3) * u(1)
      dfdu(2,4) = epsi2 * u(4) * u(1)
      dfdu(3,1) = - epsi2 * u(1) * u(4) 
      dfdu(3,2) = - epsi2 * u(2) * u(4)
      dfdu(3,3) = - 2.d0 * u(3) * u(4) 
      dfdu(3,4) = - 3.d0 * u(4) * u(4) - u(3) * u(3) - epsi
     &            * (u(1) * u(1) + u(2) * u(2))
      dfdu(4,1) = epsi2 * u(1) * u(3) 
      dfdu(4,2) = epsi2 * u(1) * u(3) 
      dfdu(4,3) = 3.d0 * u(3) * u(3) + u(4) * u(4) + epsi
     &            * (u(1) * u(1) + u(2) * u(2))
      dfdu(4,4) = 2.d0 * u(4) * u(3)
c
      dfdux(1,1) = - 0.5d0
      dfdux(1,2) = 0.d0
      dfdux(1,3) = 0.d0
      dfdux(1,4) = 0.d0
      dfdux(2,1) = 0.d0
      dfdux(2,2) = - 0.5d0
      dfdux(2,3) = 0.d0
      dfdux(2,4) = 0.d0
      dfdux(3,1) = 0.d0
      dfdux(3,2) = 0.d0
      dfdux(3,3) = 0.5d0
      dfdux(3,4) = 0.d0
      dfdux(4,1) = 0.d0
      dfdux(4,2) = 0.d0
      dfdux(4,3) = 0.d0
      dfdux(4,4) = 0.5d0
c
      dfduxx(1,1) = 0.d0
      dfduxx(1,2) = - 0.5d0
      dfduxx(1,3) = 0.d0
      dfduxx(1,4) = 0.d0
      dfduxx(2,1) = 0.5d0
      dfduxx(2,2) = 0.d0
      dfduxx(2,3) = 0.d0
      dfduxx(2,4) = 0.d0
      dfduxx(3,1) = 0.d0
      dfduxx(3,2) = 0.d0
      dfduxx(3,3) = 0.d0
      dfduxx(3,4) = - 0.5d0
      dfduxx(4,1) = 0.d0
      dfduxx(4,2) = 0.d0
      dfduxx(4,3) = 0.5d0
      dfduxx(4,4) = 0.d0
c
      return
      end
c-----------------------------------------------------------------------
      subroutine difbxa(t, u, ux, dbdu, dbdux, dbdt, npde)
c-----------------------------------------------------------------------
c Purpose:
c       The subroutine is used to define the differentiated boundary 
c       conditions at the left spatial end point x = xa. For the 
c       boundary condition equation
c                              b(t, u, ux) = 0
c       the partial derivatives db/du, db/dux, and db/dt are supplied 
c       by this routine.
c
c-----------------------------------------------------------------------
c Subroutine Parameters:
c Input:
        integer                 npde
c                               The number of PDEs in the system.
c
        double precision        t
c                               The current time coordinate.
c
        double precision        u(npde)
c                               u(1:npde) is the approximation of the
c                               solution at the point (t,x).
c
        double precision        ux(npde)
c                               ux(1:npde) is the approximation of the
c                               spatial derivative of the solution at
c                               the point (t,x).
c
c Output:
        double precision        dbdu(npde,npde)
c                               dbdu(i,j) is the partial derivative
c                               of the i-th component of the vector b
c                               with respect to the j-th component
c                               of the unknown function u.
c
        double precision        dbdux(npde,npde)
c                               dbdux(i,j) is the partial derivative
c                               of the i-th component of the vector b
c                               with respect to the j-th component
c                               of the spatial derivative of the 
c                               unknown function u.
c
        double precision        dbdt(npde)
c                               dbdt(i) is the partial derivative
c                               of the i-th component of the vector b
c                               with respect to time t.
c-----------------------------------------------------------------------
c
c     Assign dbdu(1:npde,1:npde), dbdu(1:npde,1:npde), and dbdt(1:npde)
c     according to the right boundary condition equation in terms of 
c     u(1:npde), ux(1:npde), uxx(1:npde).
c
      dbdu(1,1) = 0.d0
      dbdu(1,2) = 0.d0
      dbdu(1,3) = 0.d0
      dbdu(1,4) = 0.d0
      dbdu(2,1) = 0.d0
      dbdu(2,2) = 0.d0
      dbdu(2,3) = 0.d0
      dbdu(2,4) = 0.d0
      dbdu(3,1) = 0.d0
      dbdu(3,2) = 0.d0
      dbdu(3,3) = 0.d0
      dbdu(3,4) = 0.d0
      dbdu(4,1) = 0.d0
      dbdu(4,2) = 0.d0
      dbdu(4,3) = 0.d0
      dbdu(4,4) = 0.d0
c
      dbdux(1,1) = 1.d0
      dbdux(1,2) = 0.d0
      dbdux(1,3) = 0.d0
      dbdux(1,4) = 0.d0
      dbdux(2,1) = 0.d0
      dbdux(2,2) = 1.d0
      dbdux(2,3) = 0.d0
      dbdux(2,4) = 0.d0
      dbdux(3,1) = 0.d0
      dbdux(3,2) = 0.d0
      dbdux(3,3) = 1.d0
      dbdux(3,4) = 0.d0
      dbdux(4,1) = 0.d0
      dbdux(4,2) = 0.d0
      dbdux(4,3) = 0.d0
      dbdux(4,4) = 1.d0
c
      dbdt(1) = 0.d0
      dbdt(2) = 0.d0
      dbdt(3) = 0.d0
      dbdt(4) = 0.d0
c
      return
      end
c-----------------------------------------------------------------------
      subroutine difbxb(t, u, ux, dbdu, dbdux, dbdt, npde)
c-----------------------------------------------------------------------
c Purpose:
c       The subroutine is used to define the differentiated boundary 
c       conditions at the right spatial end point 1 = xb. For the 
c       boundary condition equation
c                              b(t, u, ux) = 0
c       the partial derivatives db/du, db/dux, and db/dt are supplied 
c       by this routine.
c
c-----------------------------------------------------------------------
c Subroutine Parameters:
c Input:
        integer                 npde
c                               The number of PDEs in the system.
c
        double precision        t
c                               The current time coordinate.
c
        double precision        u(npde)
c                               u(1:npde) is the approximation of the
c                               solution at the point (t,x).
c
        double precision        ux(npde)
c                               ux(1:npde) is the approximation of the
c                               spatial derivative of the solution at
c                               the point (t,x).
c
c Output:
        double precision        dbdu(npde,npde)
c                               dbdu(i,j) is the partial derivative
c                               of the i-th component of the vector b
c                               with respect to the j-th component
c                               of the unknown function u.
c
        double precision        dbdux(npde,npde)
c                               dbdux(i,j) is the partial derivative
c                               of the i-th component of the vector b
c                               with respect to the j-th component
c                               of the spatial derivative of the 
c                               unknown function u.
c
        double precision        dbdt(npde)
c                               dbdt(i) is the partial derivative
c                               of the i-th component of the vector b
c                               with respect to time t.
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
c     Assign dbdu(1:npde,1:npde), dbdu(1:npde,1:npde), and dbdt(1:npde)
c     according to the right boundary condition equation in terms of 
c     u(1:npde), ux(1:npde), uxx(1:npde).
c
      dbdu(1,1) = 0.d0
      dbdu(1,2) = 0.d0
      dbdu(1,3) = 0.d0
      dbdu(1,4) = 0.d0
      dbdu(2,1) = 0.d0
      dbdu(2,2) = 0.d0
      dbdu(2,3) = 0.d0
      dbdu(2,4) = 0.d0
      dbdu(3,1) = 0.d0
      dbdu(3,2) = 0.d0
      dbdu(3,3) = 0.d0
      dbdu(3,4) = 0.d0
      dbdu(4,1) = 0.d0
      dbdu(4,2) = 0.d0
      dbdu(4,3) = 0.d0
      dbdu(4,4) = 0.d0
c
      dbdux(1,1) = 1.d0
      dbdux(1,2) = 0.d0
      dbdux(1,3) = 0.d0
      dbdux(1,4) = 0.d0
      dbdux(2,1) = 0.d0
      dbdux(2,2) = 1.d0
      dbdux(2,3) = 0.d0
      dbdux(2,4) = 0.d0
      dbdux(3,1) = 0.d0
      dbdux(3,2) = 0.d0
      dbdux(3,3) = 1.d0
      dbdux(3,4) = 0.d0
      dbdux(4,1) = 0.d0
      dbdux(4,2) = 0.d0
      dbdux(4,3) = 0.d0
      dbdux(4,4) = 1.d0
c
      dbdt(1) = 0.d0
      dbdt(2) = 0.d0
      dbdt(3) = 0.d0
      dbdt(4) = 0.d0
c
      return
      end
c-----------------------------------------------------------------------
      SUBROUTINE BNDXA(T, U, UX, BVAL, NPDE)
C-----------------------------------------------------------------------
C PURPOSE:
C       THE SUBROUTINE IS USED TO DEFINE THE BOUNDARY CONDITIONS AT THE
C       LEFT SPATIAL END POINT X = XA.
C                           B(T, U, UX) = 0
C
C-----------------------------------------------------------------------
C SUBROUTINE PARAMETERS:
C INPUT:
        INTEGER                 NPDE
C                               THE NUMBER OF PDES IN THE SYSTEM.
C
        DOUBLE PRECISION        T
C                               THE CURRENT TIME COORDINATE.
C
        DOUBLE PRECISION        U(NPDE)
C                               U(1:NPDE) IS THE APPROXIMATION OF THE
C                               SOLUTION AT THE POINT (T,XA).
C
        DOUBLE PRECISION        UX(NPDE)
C                               UX(1:NPDE) IS THE APPROXIMATION OF THE
C                               SPATIAL DERIVATIVE OF THE SOLUTION AT
C                               THE POINT (T,XA).
C
C OUTPUT:
        DOUBLE PRECISION        BVAL(NPDE)
C                               BVAL(1:NPDE) IS THE BOUNDARY CONTIDITION
C                               AT THE LEFT BOUNDARY POINT.
C-----------------------------------------------------------------------
      BVAL(1) = UX(1) 
      BVAL(2) = UX(2) 
      BVAL(3) = UX(3)
      BVAL(4) = UX(4)
C
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE BNDXB(T, U, UX, BVAL, NPDE)
C-----------------------------------------------------------------------
C PURPOSE:
C       THE SUBROUTINE IS USED TO DEFINE THE BOUNDARY CONDITIONS AT THE
C       RIGHT SPATIAL END POINT X = XB.
C                           B(T, U, UX) = 0
C
C-----------------------------------------------------------------------
C SUBROUTINE PARAMETERS:
C INPUT:
        INTEGER                 NPDE
C                               THE NUMBER OF PDES IN THE SYSTEM.
C
        DOUBLE PRECISION        T
C                               THE CURRENT TIME COORDINATE.
C
        DOUBLE PRECISION        U(NPDE)
C                               U(1:NPDE) IS THE APPROXIMATION OF THE
C                               SOLUTION AT THE POINT (T,XB).
C
        DOUBLE PRECISION        UX(NPDE)
C                               UX(1:NPDE) IS THE APPROXIMATION OF THE
C                               SPATIAL DERIVATIVE OF THE SOLUTION AT
C                               THE POINT (T,XB).
C
C OUTPUT:
        DOUBLE PRECISION        BVAL(NPDE)
C                               BVAL(1:NPDE) IS THE BOUNDARY CONTIDITION
C                               AT THE RIGHT BOUNDARY POINT.
C-----------------------------------------------------------------------
      BVAL(1) = UX(1)
      BVAL(2) = UX(2)
      BVAL(3) = UX(3)
      BVAL(4) = UX(4)
C
      RETURN
      END
C-----------------------------------------------------------------------
      subroutine f(t, x, u, ux, uxx, fval, npde)
c-----------------------------------------------------------------------
c Purpose:
c       This subroutine defines the right hand side vector of the 
c       npde dimensional parabolic partial differential equation 
c                        ut = f(t, x, u, ux, uxx). 
c
c-----------------------------------------------------------------------
c Subroutine Parameters:
c Input:
        integer                 npde
c                               The number of PDEs in the system.
c
        double precision        t
c                               The current time coordinate.
c
        double precision        x
c                               The current spatial coordinate.
c
        double precision        u(npde)
c                               u(1:npde) is the approximation of the
c                               solution at the point (t,x).
c
        double precision        ux(npde)
c                               ux(1:npde) is the approximation of the
c                               spatial derivative of the solution at
c                               the point (t,x).
c
        double precision        uxx(npde)
c                               uxx(1:npde) is the approximation of the
c                               second spatial derivative of the 
c                               solution at the point (t,x).
c
c Output:
        double precision        fval(npde)
c                               fval(1:npde) is the right hand side
c                               vector f(t, x, u, ux, uxx) of the PDE.
c-----------------------------------------------------------------------
c
c     Assign fval(1:npde) according to the right hand side of the PDE 
c     in terms of u(1:npde), ux(1:npde), uxx(1:npde).
c
      fval(1) = - 0.5d0 * ux(1) - 0.5d0 * uxx(2) - u(2)
     &          * ((u(1) * u(1) + u(2) * u(2)) + 2.d0/3.d0
     &          * ((u(3) * u(3) + u(4) * u(4))))
      fval(2) = - 0.5d0 * ux(2) + 0.5d0 * uxx(1) + u(1)
     &          * ((u(1) * u(1) + u(2) * u(2)) + 2.d0/3.d0
     &          * ((u(3) * u(3) + u(4) * u(4))))
      fval(3) = 0.5d0 * ux(3) - 0.5d0 * uxx(4) - u(4)
     &          * ((u(3) * u(3) + u(4) * u(4)) + 2.d0/3.d0
     &          * ((u(1) * u(1) + u(2) * u(2))))
      fval(4) = 0.5d0 * ux(4) + 0.5d0 * uxx(3) + u(3)
     &          * ((u(3) * u(3) + u(4) * u(4)) + 2.d0/3.d0
     &          * ((u(1) * u(1) + u(2) * u(2))))
c
      return
      end
c-----------------------------------------------------------------------
      subroutine uinit(x, u, npde)
c-----------------------------------------------------------------------
c Purpose:
c       This subroutine is used to return the npde-vector of initial 
c       conditions of the unknown function at the initial time t = t0 
c       at the spatial coordinate x.
c
c-----------------------------------------------------------------------
c Subroutine Parameters:
c Input:
        double precision        x
c                               The spatial coordinate.
c
        integer                 npde
c                               The number of PDEs in the system.
c
c Output:
        double precision        u(npde)
c                               u(1:npde) is vector of initial values of
c                               the unknown function at t = t0 and the
c                               given value of x.
c-----------------------------------------------------------------------
        double precision        tempt1, tempt2

c
c     Assign u(1:npde) the initial values of u(t0,x).
c
      
      tempt1 = sqrt(6.d0/5.d0)
      tempt2 = sqrt(2.d0)

      u(1) = tempt1 / cosh(tempt2 * x) * cos(0.5d0 * x)
      u(2) = tempt1 / cosh(tempt2 * x) * sin(0.5d0 * x)  
      u(3) = tempt1 / cosh(tempt2 * x) * cos(1.5d0 * x)  
      u(4) = tempt1 / cosh(tempt2 * x) * sin(1.5d0 * x)  
c
      return
      end

