```!   D03PFF Example Program Text
!   Mark 25 Release. NAG Copyright 2014.

Module d03pffe_mod

!     D03PFF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Accessibility Statements ..
Private
Public                               :: bndry1, bndry2, exact, nmflx1,   &
nmflx2, pdedef
!     .. Parameters ..
Integer, Parameter, Public           :: nin = 5, nout = 6, npde1 = 2,    &
npde2 = 1
Contains
Subroutine exact(t,u,npde,x,npts)
!       Exact solution (for comparison and b.c. purposes)

!       .. Use Statements ..
Use nag_library, Only: x01aaf
!       .. Parameters ..
Real (Kind=nag_wp), Parameter        :: half = 0.5_nag_wp
Real (Kind=nag_wp), Parameter        :: two = 2.0_nag_wp
!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In)      :: t
Integer, Intent (In)                 :: npde, npts
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out)     :: u(npde,npts)
Real (Kind=nag_wp), Intent (In)      :: x(npts)
!       .. Local Scalars ..
Real (Kind=nag_wp)                   :: pi, px1, px2, x1, x2
Integer                              :: i
!       .. Intrinsic Procedures ..
Intrinsic                            :: exp, sin
!       .. Executable Statements ..
pi = x01aaf(pi)

Do i = 1, npts
x1 = x(i) + t
x2 = x(i) - 3.0_nag_wp*t
px1 = half*sin(two*pi*x1**2)
px2 = half*sin(two*pi*x2**2)
u(1,i) = half*(exp(x2)+exp(x1)+px2-px1) - t*(x1+x2)
u(2,i) = (exp(x2)-exp(x1)+px2+px1) + x1*x2 + 8.0_nag_wp*t**2
End Do
Return
End Subroutine exact
Subroutine bndry1(npde,npts,t,x,u,ibnd,g,ires)

!       .. Parameters ..
Real (Kind=nag_wp), Parameter        :: one = 1.0_nag_wp
Real (Kind=nag_wp), Parameter        :: two = 2.0_nag_wp
!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In)      :: t
Integer, Intent (In)                 :: ibnd, npde, npts
Integer, Intent (Inout)              :: ires
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out)     :: g(npde)
Real (Kind=nag_wp), Intent (In)      :: u(npde,3), x(npts)
!       .. Local Scalars ..
Real (Kind=nag_wp)                   :: c, exu1, exu2
!       .. Local Arrays ..
Real (Kind=nag_wp)                   :: ue(2,1)
!       .. Executable Statements ..
If (ibnd==0) Then
Call exact(t,ue,npde,x(1),1)
c = (x(2)-x(1))/(x(3)-x(2))
exu1 = (one+c)*u(1,2) - c*u(1,3)
exu2 = (one+c)*u(2,2) - c*u(2,3)
g(1) = two*u(1,1) + u(2,1) - two*ue(1,1) - ue(2,1)
g(2) = two*u(1,1) - u(2,1) - two*exu1 + exu2
Else
Call exact(t,ue,npde,x(npts),1)
c = (x(npts)-x(npts-1))/(x(npts-1)-x(npts-2))
exu1 = (one+c)*u(1,2) - c*u(1,3)
exu2 = (one+c)*u(2,2) - c*u(2,3)
g(1) = two*u(1,1) - u(2,1) - two*ue(1,1) + ue(2,1)
g(2) = two*u(1,1) + u(2,1) - two*exu1 - exu2
End If
Return
End Subroutine bndry1
Subroutine nmflx1(npde,t,x,uleft,uright,flux,ires)

!       .. Parameters ..
Real (Kind=nag_wp), Parameter        :: half = 0.5_nag_wp
!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In)      :: t, x
Integer, Intent (Inout)              :: ires
Integer, Intent (In)                 :: npde
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out)     :: flux(npde)
Real (Kind=nag_wp), Intent (In)      :: uleft(npde), uright(npde)
!       .. Local Scalars ..
Real (Kind=nag_wp)                   :: ltmp, rtmp
!       .. Executable Statements ..
ltmp = 3.0_nag_wp*uleft(1) + 1.5_nag_wp*uleft(2)
rtmp = uright(1) - half*uright(2)
flux(1) = half*(ltmp-rtmp)
flux(2) = ltmp + rtmp
Return
End Subroutine nmflx1
Subroutine pdedef(npde,t,x,u,ux,p,c,d,s,ires)

!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In)      :: t, x
Integer, Intent (Inout)              :: ires
Integer, Intent (In)                 :: npde
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out)     :: c(npde), d(npde),              &
p(npde,npde), s(npde)
Real (Kind=nag_wp), Intent (In)      :: u(npde), ux(npde)
!       .. Executable Statements ..
p(1,1) = 1.0_nag_wp
c(1) = 0.01_nag_wp
d(1) = ux(1)
s(1) = u(1)
Return
End Subroutine pdedef
Subroutine bndry2(npde,npts,t,x,u,ibnd,g,ires)

!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In)      :: t
Integer, Intent (In)                 :: ibnd, npde, npts
Integer, Intent (Inout)              :: ires
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out)     :: g(npde)
Real (Kind=nag_wp), Intent (In)      :: u(npde,3), x(npts)
!       .. Executable Statements ..
If (ibnd==0) Then
g(1) = u(1,1) - 3.0_nag_wp
Else
g(1) = u(1,1) - 5.0_nag_wp
End If
Return
End Subroutine bndry2
Subroutine nmflx2(npde,t,x,uleft,uright,flux,ires)

!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In)      :: t, x
Integer, Intent (Inout)              :: ires
Integer, Intent (In)                 :: npde
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out)     :: flux(npde)
Real (Kind=nag_wp), Intent (In)      :: uleft(npde), uright(npde)
!       .. Executable Statements ..
If (x>=0.0E0_nag_wp) Then
flux(1) = x*uleft(1)
Else
flux(1) = x*uright(1)
End If
Return
End Subroutine nmflx2
End Module d03pffe_mod
Program d03pffe

!     D03PFF Example Main Program

!     .. Use Statements ..
Use d03pffe_mod, Only: nout
!     .. Implicit None Statement ..
Implicit None
!     .. Executable Statements ..
Write (nout,*) 'D03PFF Example Program Results'

Call ex1

Call ex2

Contains
Subroutine ex1

!       .. Use Statements ..
Use nag_library, Only: d03pff, d03pfp, nag_wp
Use d03pffe_mod, Only: bndry1, exact, nin, nmflx1, npde1
!       .. Local Scalars ..
Real (Kind=nag_wp)                   :: dx, tout, ts, tsmax
Integer                              :: i, ifail, inc, ind, it, itask, &
itrace, lisave, lrsave,        &
nfuncs, niters, njacs, nop,    &
npde, npts, nsteps, outpts
!       .. Local Arrays ..
Real (Kind=nag_wp)                   :: acc(2)
Real (Kind=nag_wp), Allocatable      :: rsave(:), u(:,:), ue(:,:),     &
x(:), xout(:)
Integer, Allocatable                 :: isave(:)
!       .. Intrinsic Procedures ..
Intrinsic                            :: real
!       .. Executable Statements ..
Write (nout,*)
Write (nout,*)
Write (nout,*) 'Example 1'
Write (nout,*)
!       Skip heading in data file
npde = npde1
Read (nin,*) npts, inc, outpts
lisave = 24 + npde*npts
lrsave = (11+9*npde)*npde*npts + (32+3*npde)*npde + 7*npts + 54

Allocate (rsave(lrsave),u(npde,npts),ue(npde,outpts),x(npts), &
xout(outpts),isave(lisave))

!       Initialise mesh
dx = 1.0_nag_wp/real(npts-1,kind=nag_wp)
Do i = 1, npts
x(i) = real(i-1,kind=nag_wp)*dx
End Do

!       Set initial values
Call exact(ts,u,npde,x,npts)

ind = 0

Write (nout,99992) npts
Write (nout,99991) acc(1)
Write (nout,99990) acc(2)
Write (nout,99999)

Do it = 1, 2
tout = 0.1E0_nag_wp*real(it,kind=nag_wp)

!         ifail: behaviour on error exit
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d03pff(npde,ts,tout,d03pfp,nmflx1,bndry1,u,npts,x,acc,tsmax, &

!         Set output points

nop = 0
Do i = 1, npts, inc
nop = nop + 1
xout(nop) = x(i)
End Do

!         Check against exact solution

Call exact(tout,ue,npde,xout,nop)
nop = 1
i = 1
Write (nout,99998) ts, xout(nop), u(1,i), ue(1,nop), u(2,i), &
ue(2,nop)
Do i = inc + 1, npts, inc
nop = nop + 1
Write (nout,99997) xout(nop), u(1,i), ue(1,nop), u(2,i), ue(2,nop)
End Do
End Do

!       Print integration statistics (reasonably rounded)
nsteps = 5*((isave(1)+2)/5)
nfuncs = 50*((isave(2)+25)/50)
njacs = 10*((isave(3)+5)/10)
niters = 50*((isave(5)+25)/50)
Write (nout,99996) nsteps
Write (nout,99995) nfuncs
Write (nout,99994) njacs
Write (nout,99993) niters

Return

99999   Format (/5X,'t',9X,'x',8X,'Approx u',4X,'Exact u',5X,'Approx v',4X, &
'Exact v')
99998   Format (/1X,F6.3,5(3X,F9.4))
99997   Format (7X,5(3X,F9.4))
99996   Format (/' Number of time steps           (nearest 5)  = ',I6)
99995   Format (' Number of function evaluations (nearest 50) = ',I6)
99994   Format (' Number of Jacobian evaluations (neaerst 10) = ',I6)
99993   Format (' Number of iterations           (nearest 50) = ',I6)
99992   Format (/' Number of mesh points used = ',I4)
99991   Format (' Relative tolerance used    = ',E10.3)
99990   Format (' Absolute tolerance used    = ',E10.3)
End Subroutine ex1
Subroutine ex2

!       .. Use Statements ..
Use nag_library, Only: d03pff, nag_wp
Use d03pffe_mod, Only: bndry2, nin, nmflx2, npde2, pdedef
!       .. Local Scalars ..
Real (Kind=nag_wp)                   :: dx, tout, ts, tsmax
Integer                              :: i, ifail, ind, it, itask,      &
itrace, lisave, lrsave,        &
nfuncs, niters, njacs, npde,   &
npts, nsteps, outpts
!       .. Local Arrays ..
Real (Kind=nag_wp)                   :: acc(2)
Real (Kind=nag_wp), Allocatable      :: rsave(:), u(:,:), x(:), xout(:)
Integer, Allocatable                 :: iout(:), isave(:)
!       .. Intrinsic Procedures ..
Intrinsic                            :: real
!       .. Executable Statements ..
Write (nout,*)
Write (nout,*)
Write (nout,*) 'Example 2'
Write (nout,*)
npde = npde2
Read (nin,*) npts, outpts
lisave = 24 + npde*npts
lrsave = (11+9*npde)*npde*npts + (32+3*npde)*npde + 7*npts + 54

Allocate (rsave(lrsave),u(npde,npts),x(npts),xout(outpts), &
isave(lisave),iout(outpts))

!       Initialise mesh
dx = 2.0_nag_wp/real(npts-1,kind=nag_wp)
Do i = 1, npts
x(i) = -1.0_nag_wp + real(i-1,kind=nag_wp)*dx
End Do

!       Set initial values
u(1,1:npts) = x(1:npts) + 4.0_nag_wp

ind = 0

!       Set output points from inout indices
Do i = 1, outpts
xout(i) = x(iout(i))
End Do

!       Two output value of t: tout read from file and 10.0

Read (nin,*) ts, tout

Write (nout,99995) npts
Write (nout,99994) acc(1)
Write (nout,99993) acc(2)
Write (nout,99992) xout(1:outpts)

Do it = 1, 2

ifail = 0
Call d03pff(npde,ts,tout,pdedef,nmflx2,bndry2,u,npts,x,acc,tsmax, &

Write (nout,99991) ts, (u(1,iout(i)),i=1,outpts)
tout = 10.0_nag_wp

End Do

!       Print integration statistics (reasonably rounded)
nsteps = 5*((isave(1)+2)/5)
nfuncs = 50*((isave(2)+25)/50)
njacs = 10*((isave(3)+5)/10)
niters = 50*((isave(5)+25)/50)
Write (nout,99999) nsteps
Write (nout,99998) nfuncs
Write (nout,99997) njacs
Write (nout,99996) niters

Return

99999   Format (/' Number of time steps           (nearest 5)  = ',I6)
99998   Format (' Number of function evaluations (nearest 50) = ',I6)
99997   Format (' Number of Jacobian evaluations (neaerst 10) = ',I6)
99996   Format (' Number of iterations           (nearest 50) = ',I6)
99995   Format (/' Number of mesh points used = ',I4)
99994   Format (' Relative tolerance used    = ',E10.3)
99993   Format (' Absolute tolerance used    = ',E10.3)
99992   Format (/1X,'x',10X,7F9.4)
99991   Format (1X,'u(t=',F6.3,')',7F9.4)
End Subroutine ex2
End Program d03pffe
```