C     EEE2022-Ex5

! U(1) = u(x,t)
! U(2) = w(x,t) [Backward Integral from x to 1]
! U(3) = v(x,t) [Forward Integral from 0 to x]

      SUBROUTINE F(T, X, U, UX, UXX, FVAL, NPDE)
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
      PI = 4.0D0 * DATAN(1.0D0)

C     Heat Equation
      FVAL(1) = UXX(1) + (PI**2 - 1.0D0) * EXP(-T) *
     #          (DSIN(PI*X) + DCOS(PI*X))

C     ODE for w: w_x = -sin(pi*x)*u in residual form
      FVAL(2) = UX(2) + DSIN(PI*X) * U(1)

C     ODE for v: v_x = cos(pi*x)*u in residual form
      FVAL(3) = UX(3) - DCOS(PI*X) * U(1)
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-----------------------------------------------------------------------
C
C     PDE boundary condition in residual form: u - 2 * w = 0
      BVAL(1) = U(1) - 2.0D0 * U(2)

C     Left boundary condition for v: v = 0
      BVAL(2) = U(3)

C     There is no left BC for w so use
C     Residual of ODE for w, evaluated at x=0: w_x + sin(pi*0)*u = 0 => w_x = 0
      BVAL(3) = UX(2)
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-----------------------------------------------------------------------
C
      PI = 4.0D0 * DATAN(1.0D0)
C
C     Residual form of left BC of PDE: u = -2v = 0
      BVAL(1) = U(1) + 2.0D0 * U(3)

C     Right BC for w: w = 0
      BVAL(2) = U(2)

C     There is no right BC for v so use the ODE itself, evaluated at x=1:
C     Residual v_x - cos(pi*1)*u = 0 => v_x + u = 0
      BVAL(3) = UX(3) + U(1)
C
      RETURN
      END

      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
      PI = 4.0D0 * DATAN(1.0D0)
C
C     Initial Condition for u: u(x,0) = sin(pi*x) + cos(pi*x)
      U(1) = DSIN(PI*X) + DCOS(PI*X)

C     Initial Condition for w: w(x,0) = 0.5 - (x/2 - sin(2*pi*x)/(4*pi) + sin(pi*x)**2/(2*pi))
      U(2) = 0.5D0 - (X/2.0D0 - DSIN(2.0D0*PI*X)/(4.0D0*PI) + 
     #       (DSIN(PI*X)**2)/(2.0D0*PI))
C
C     Initial Condition for v: v(x,0) = sin(pi*x)**2/(2*pi) + x/2 + sin(2*pi*x)/(4*pi)
      U(3) = (DSIN(PI*X)**2)/(2.0D0*PI) + X/2.0D0 + 
     #       DSIN(2.0D0*PI*X)/(4.0D0*PI)
C
      RETURN
      END

      SUBROUTINE TRUU(T, X, U, NPDE)
C-----------------------------------------------------------------------
C PURPOSE:
C     THIS FUNCTION PROVIDES THE EXACT SOLUTIONS
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
C OUTPUT:
        DOUBLE PRECISION        U(NPDE)
C                               U(1:NPDE) IS THE EXACT SOLUTION AT THE
C                               POINT (T,X).
C-----------------------------------------------------------------------

      DOUBLE PRECISION PI, ET
C
      PI = 4.0D0 * DATAN(1.0D0)
      ET = DEXP(-T)

C     Exact solution for u
      U(1) = ET * (DSIN(PI*X) + DCOS(PI*X))

C     Exact solution for w
      U(2) = ET * (0.5D0 - (X/2.0D0 - DSIN(2.0D0*PI*X)/(4.0D0*PI) + 
     #       (DSIN(PI*X)**2)/(2.0D0*PI)))

C     Exact solution for v
      U(3) = ET * ((DSIN(PI*X)**2)/(2.0D0*PI) + X/2.0D0 + 
     #       DSIN(2.0D0*PI*X)/(4.0D0*PI))
C
      RETURN
      END