C     AlYo2023-Ex1

C     U(1) = u, U(2) = v, U(3) = w
C-----------------------------------------------------------------------
      SUBROUTINE F(T, X, U, UX, UXX, FVAL, NPDE)
C-----------------------------------------------------------------------
        INTEGER          NPDE
        DOUBLE PRECISION T, X, U(NPDE), UX(NPDE), UXX(NPDE), FVAL(NPDE)

        ! PDE for u
        FVAL(1) = UXX(1) - U(1)

        ! ODE for v
        FVAL(2) = UX(2) - X * U(1)

        ! ODE for w
        FVAL(3) = UX(3) - (1.0D0 - X) * U(1)

        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-----------------------------------------------------------------------
        INTEGER                 NPDE
        DOUBLE PRECISION        T, U(NPDE), UX(NPDE), BVAL(NPDE)

        ! Residual form of the ODE for w at x=0
        BVAL(1) = UX(3) -  U(1)

        ! Boundary condition for v
        BVAL(2) = U(2)

        ! Boundary condition for w
        BVAL(3) = U(3)

        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-----------------------------------------------------------------------
        INTEGER                 NPDE
        DOUBLE PRECISION        T, U(NPDE), UX(NPDE), BVAL(NPDE)
        REAL(8), PARAMETER :: PI = 3.141592653589793D0

        ! Reformuated integral constraint involving v
        BVAL(1) = U(2) - PI * EXP(-2.0D0 * T)

        ! Reformuated integral constraint involving w
        BVAL(2) = U(3) - (2.0D0 - PI) * EXP(-2.0D0 * T)

        ! Residual form of the ODE for v at x=pi
        BVAL(3) = UX(2) - PI * U(1)

        RETURN
        END

C-----------------------------------------------------------------------
      SUBROUTINE UINIT(X, U, NPDE)
C-----------------------------------------------------------------------
        INTEGER                 NPDE
        DOUBLE PRECISION        X, U(NPDE)

        ! Initial condition for u
        U(1) = SIN(X)

        ! Initial condition for v
        U(2) = SIN(X) - X * COS(X)

        ! Initial condition for w
        U(3) = (1.0D0 - X) * (-COS(X)) - SIN(X) + 1.0D0

        RETURN
        END

C-----------------------------------------------------------------------
      SUBROUTINE TRUU(T, X, U, NPDE)
C-----------------------------------------------------------------------
        INTEGER                 NPDE
        DOUBLE PRECISION        T, X, U(NPDE)
        REAL(8) :: ET

        ET = EXP(-2.0D0 * T)

        ! Exact solution for u
        U(1) = ET * SIN(X)

        ! Exact solution for v
        U(2) = ET * (SIN(X) - X * COS(X))

        ! Exact solution for w
        U(3) = ET * ((1.0D0 - X) * (-COS(X)) - SIN(X) + 1.0D0)

        RETURN
      END