c     ZoSh2013-Ex1

! 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]

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

C     Heat Equation
      FVAL(1) = UXX(1) 

C     ODE for w: w_x = -(x + t)*u in residual form
      FVAL(2) = UX(2) + (X + T) * U(1)

C     ODE for v: v_x = (t * exp(-x))*u in residual form
      FVAL(3) = UX(3) - (T * EXP(-X)) * U(1)

      RETURN
      END


      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-----------------------------------------------------------------------

      REAL(8) :: RHS1
   
      RHS1 = 0.5D0 + EXP(-T) - T - 
     #      EXP(-T) * (-1.0D0 + COS(1.0D0) + SIN(1.0D0) + T*SIN(1.0D0))

C       Residual form of left BC of PDE: u = w + RHS1
        BVAL(1) = U(1) - U(2) - RHS1

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

C       There is no left BC for w so use the ODE itself:
C       Residual form of ODE: w_x + (x+t)u = 0, evaluated at x=0
        BVAL(3) = UX(2) + T * U(1)

      RETURN
      END

      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-----------------------------------------------------------------------

      REAL(8) :: RHS2

      RHS2 = 1.0D0 + EXP(-T)*COS(1.0D0) - 
     #      (T / (2.0D0*EXP(1.0D0))) * 
     #     (2.0D0*(EXP(1.0D0) - 1.0D0) + EXP(-T)*(EXP(1.0D0) 
     #     - COS(1.0D0) + SIN(1.0D0)))

C       Residual form of right BC of PDE: u = v + RHS2
        BVAL(1) = U(1) - U(3) - RHS2

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

C       There is no right BC for v so use the ODE itself:
C       Residual form of v_x - (t*exp(-x))u = 0, evaluated at x=1
        BVAL(3) = UX(3) - (T / EXP(1.0D0)) * U(1)

      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.
        REAL(8), PARAMETER :: PI = 4.d0*datan(1.d0)

c-----------------------------------------------------------------------

C
C     ASSIGN U(1:NPDE) THE INITIAL VALUES OF U(T0,X).
C

    
C     Initial Condition for u
      U(1) = 1.0D0 + COS(X)
    
C     Initial Condition for w 
      U(2) = (0.5D0 + SIN(1.0D0) + COS(1.0D0)) - 
     #      (0.5D0*X**2 + X*SIN(X) + COS(X))
           
C     Initial Condition for v 
C     Since this kernel is (t * exp(-s)), at t=0, this is 0
      U(3) = 0.0D0

      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-----------------------------------------------------------------------

      REAL(8) :: ET, TET

      ET = EXP(-T)
      TET = T * ET

C     Exact solution for u
      U(1) = 1.0D0 + ET * COS(X)

C     Exact solution for w
      U(2) = (0.5D0*(1.0D0 - X**2) + T*(1.0D0 - X)) + 
     #      ET * ((SIN(1.0D0) + COS(1.0D0)) - (X*SIN(X) + COS(X))) + 
     #      TET * (SIN(1.0D0) - SIN(X))

C     Exact solution for v
      U(3) = T * (1.0D0 - EXP(-X)) + 
     #      0.5D0 * TET * (EXP(-X)*(SIN(X) - COS(X)) + 1.0D0)

      RETURN
      END
