! pde45h_eq.f90  homogeneous Biharmonic equation in five dimensions
! solve Uwwww(w,x,y,z,t) + Uxxxx(w,x,y,z,t) + Uyyyy(w,x,y,z,t) + 
!       Uzzzz(w,x,y,z,t) + Utttt(w,x,y,z,t) + 2 Uww(w,x,y,z,t) + 
!       2 Uxx(w,x,y,z,t) + 2 Uyy(w,x,y,z,t) + 2 Uzz(w,x,y,z,t) +
!       2 Utt(w,x,y,z,t) + 5 U(w,x,y,z,t) = 0
!
! boundary conditions computed using U(w,x,y,z,t)
! analytic solution is U(w,x,y,z,t) = sin(w+x+y+z+t)
!

module global
  implicit none
  ! problem parameters
  integer, parameter :: nw=6
  integer, parameter :: nx=6
  integer, parameter :: ny=6
  integer, parameter :: nz=6
  integer, parameter :: nt=6
  integer, parameter :: nwxyzt=(nw-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2)
  double precision :: wmin = 0.0
  double precision :: wmax = 1.5708
  double precision :: xmin = 0.0
  double precision :: xmax = 1.5708
  double precision :: ymin = 0.0
  double precision :: ymax = 1.5708
  double precision :: zmin = 0.0
  double precision :: zmax = 1.5708
  double precision :: tmin = 0.0
  double precision :: tmax = 1.5708
  double precision :: hw, hx, hy, hz, ht
  double precision, dimension(nwxyzt,nwxyzt) :: A ! stiffness matrix
     ! A(s(i,ii,iii,iiii,iiiii),s(j,jj,jjj,jjjj,jjjjj))
  double precision, dimension(nw) :: wg ! w grid coordinates
  double precision, dimension(nx) :: xg ! x
  double precision, dimension(ny) :: yg ! y
  double precision, dimension(nz) :: zg ! z
  double precision, dimension(nt) :: tg ! t
  double precision, dimension(nwxyzt) :: R ! RHS
  double precision, dimension(nwxyzt) :: U ! unknown solution
  double precision, dimension(nwxyzt) :: Ua ! known solution for checking

  interface
    function uana(w,x,y,z,t) result(v) ! analytic solution and boundary
      double precision, intent(in) :: w,x,y,z,t
      double precision :: v
    end function uana
    subroutine simeq(n, sz, A, Y, X)
      integer, intent(in) :: n  ! number of equations
      integer, intent(in) :: sz ! dimension of arrays
      double precision, dimension(sz,sz), intent(in) :: A
      double precision, dimension(sz), intent(in) :: Y
      double precision, dimension(sz), intent(out) :: X
    end subroutine simeq
    subroutine rderiv1(order, npoints, point, hx, cx)
      integer, intent(in) :: order, npoints, point
      double precision, intent(in) :: hx
      double precision, dimension(npoints), intent(out) :: cx
    end subroutine rderiv1
    function s(i, ii, iii, iiii, iiiii) result(v)
      integer, intent(in) :: i, ii, iii, iiii, iiiii
      integer :: v
    end function s
    subroutine buildA()
    end subroutine buildA
    function eval_derivative(word, xord, yord, zord, tord, &
                             i, ii, iii, iiii, iiiii, UU) result (discrete)
      integer, intent(in) :: word, xord, yord, zord, tord
      integer, intent(in) :: i, ii, iii, iiii, iiiii
      double precision, dimension(*), intent(in) :: UU
      double precision :: discrete
    end function eval_derivative
    subroutine check_soln(UU) ! for specific PDE
      double precision, dimension(*), intent(in) :: UU
    end subroutine check_soln
    subroutine write_soln(UU) ! for specific PDE
      double precision, dimension(*), intent(in) :: UU
    end subroutine write_soln

  end interface

end module global

! functions called by program
  function uana(w,x,y,z,t) result(v) ! analytic solution and boundary
    implicit none
    double precision, intent(in) :: w,x,y,z,t
    double precision :: v

    v =  sin(w+x+y+z+t)
  end function uana

  function s(i, ii, iii, iiii, iiiii) result(v)
    use global
    implicit none
    integer, intent(in) :: i, ii, iii, iiii, iiiii
    integer :: v

    v = (i-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2) + &
        (ii-2)*(ny-2)*(nz-2)*(nt-2) + (iii-2)*(nz-2)*(nt-2) + &
        (iiii-2)*(nt-2) + iiiii-1
  end function s

  subroutine buildA() ! dpends on specific PDE
    use global
    implicit none
    integer :: i, j, ii, jj, iii, jjj, iiii, jjjj, iiiii, jjjjj, k
    double precision :: val = 0
    double precision, dimension(nw) :: cww, cwwww
    double precision, dimension(nx) :: cxx, cxxxx
    double precision, dimension(ny) :: cyy, cyyyy
    double precision, dimension(nz) :: czz, czzzz
    double precision, dimension(nt) :: ctt, ctttt

    ! compute entries=A matrix, inside boundary
    print *, 'compute A matrix and Y RHS'


    print *, 'check rderiv1(4,nw,1,hw,cwwww)'
    call rderiv1(4,nw,1,hw,cwwww)
    do i=1,nw 
      print *, 'cwwww(',i,')=',cwwww(i)
    end do
    print *, 'check rderiv1(4,nw,nw,hw,cwwww)'
    call rderiv1(4,nw,nw,hw,cwwww)
    do i=1,nw
      print *, 'cwwww(',i,')=',cwwww(i)
    end do
    print *, ' '

    val = 0.0
    do i=2,nw-1
      do ii=2,nx-1
        do iii=2,ny-1
          do iiii=2,nz-1
            do iiiii=2,nt-1
              k = s(i,ii,iii,iiii,iiiii) ! row index

              ! do each term

              ! do d^4u/dw^4
              call rderiv1(4,nw,i,hw,cwwww)
              do j=1,nw
                val = cwwww(j)
                if(j.eq.1 .or. j.eq.nw) then
                  R(k) = R(k)-val*uana(wg(j),xg(ii),yg(iii),zg(iiii),tg(iiiii))
                else
                  A(k,s(j,ii,iii,iiii,iiiii))=A(k,s(j,ii,iii,iiii,iiiii)) + val
                end if
              end do ! j

              ! do d^4u/dx^4
              call rderiv1(4,nx,ii,hx,cxxxx)
              do jj=1,nx
                val = cxxxx(jj)
                if(jj.eq.1 .or. jj.eq.nx) then
                  R(k) = R(k)-val*uana(wg(i),xg(jj),yg(iii),zg(iiii),tg(iiiii))
                else
                  A(k,s(i,jj,iii,iiii,iiiii))=A(k,s(i,jj,iii,iiii,iiiii)) + val
                end if
              end do ! jj

              ! do d^4u/dy^4
              call rderiv1(4,ny,iii,hy,cyyyy)
              do jjj=1,ny
                val = cyyyy(jjj)
                if(jjj.eq.1 .or. jjj.eq.ny) then
                  R(k) = R(k)-val*uana(wg(i),xg(ii),yg(jjj),zg(iiii),tg(iiiii))
                else
                  A(k,s(i,ii,jjj,iiii,iiiii))=A(k,s(i,ii,jjj,iiii,iiiii)) + val
                end if
              end do ! jjj

              ! do d^4u/dz^4
              call rderiv1(4,nz,iiii,hz,czzzz)
              do jjjj=1,nz
                val = czzzz(jjjj)
                if(jjjj.eq.1 .or. jjjj.eq.nz) then
                  R(k) = R(k)-val*uana(wg(i),xg(ii),yg(iii),zg(jjjj),tg(iiiii))
                else
                  A(k,s(i,ii,iii,jjjj,iiiii))=A(k,s(i,ii,iii,jjjj,iiiii)) + val
                end if
              end do ! jjjj

              ! do d^4u/dt^4
              call rderiv1(4,nt,iiiii,ht,ctttt)
              do jjjjj=1,nt
                val = ctttt(jjjjj)
                if(jjjjj.eq.1 .or. jjjjj.eq.nt) then
                  R(k) = R(k)-val*uana(wg(i),xg(ii),yg(iii),zg(iiii),tg(jjjjj))
                else
                  A(k,s(i,ii,iii,iiii,jjjjj))=A(k,s(i,ii,iii,iiii,jjjjj)) + val
                end if
              end do ! jjjjj

              ! do 2 d^2u/dw^2
              call rderiv1(2,nw,i,hw,cww)
              do j=1,nw
                val = 2.0*cww(j)
                if(j.eq.1 .or. j.eq.nw) then
                  R(k) = R(k)-val*uana(wg(j),xg(ii),yg(iii),zg(iiii),tg(iiiii))
                else
                  A(k,s(j,ii,iii,iiii,iiiii))=A(k,s(j,ii,iii,iiii,iiiii)) + val
                end if
              end do ! j

              ! do 2 d^2u/dx^2
              call rderiv1(2,nx,ii,hx,cxx)
              do jj=1,nx
                val = 2.0*cxx(jj)
                if(jj.eq.1 .or. jj.eq.nx) then
                  R(k) = R(k)-val*uana(wg(i),xg(jj),yg(iii),zg(iiii),tg(iiiii))
                else
                  A(k,s(i,jj,iii,iiii,iiiii))=A(k,s(i,jj,iii,iiii,iiiii)) + val
                end if
              end do ! jj

              ! do 2 d^2u/dy^2
              call rderiv1(2,ny,iii,hy,cyy)
              do jjj=1,ny
                val = 2.0*cyy(jjj)
                if(jjj.eq.1 .or. jjj.eq.ny) then
                  R(k) = R(k)-val*uana(wg(i),xg(ii),yg(jjj),zg(iiii),tg(iiiii))
                else
                  A(k,s(i,ii,jjj,iiii,iiiii))=A(k,s(i,ii,jjj,iiii,iiiii)) + val
                end if
              end do ! jjj

              ! do 2 d^2u/dz^2
              call rderiv1(2,nz,iiii,hz,czz)
              do jjjj=1,nz
                val = 2.0*czz(jjjj)
                if(jjjj.eq.1 .or. jjjj.eq.nz) then
                  R(k) = R(k)-val*uana(wg(i),xg(ii),yg(iii),zg(jjjj),tg(iiiii))
                else
                  A(k,s(i,ii,iii,jjjj,iiiii))=A(k,s(i,ii,iii,jjjj,iiiii)) + val
                end if
              end do ! jjjj

              ! do 2 d^2u/dt^2
              call rderiv1(2,nt,iiiii,ht,ctt)
              do jjjjj=1,nt
                val = 2.0*ctt(jjjjj)
                if(jjjjj.eq.1 .or. jjjjj.eq.nt) then
                  R(k) = R(k)-val*uana(wg(i),xg(ii),yg(iii),zg(iiii),tg(jjjjj))
                else
                  A(k,s(i,ii,iii,iiii,jjjjj))=A(k,s(i,ii,iii,iiii,jjjjj)) + val
                end if
              end do ! jjjjj

              ! do 5 u(w,x,y,z,t)
              A(k,s(i,ii,iii,iiii,iiiii))=A(k,s(i,ii,iii,iiii,iiiii)) + 5.0

              ! do f(w,x,y,z,t) RHS, none
              ! R(k) = R(k) + f(w,x,y,z,t)

              ! end terms

            end do ! iiiii
          end do ! iiii
        end do ! iii
      end do ! ii 
    end do ! i
  end subroutine buildA

  function eval_derivative(word, xord, yord, zord, tord, &
                           i, ii, iii, iiii, iiiii, UU) result (discrete)
    use global
    implicit none
    integer, intent(in) :: word, xord, yord, zord, tord
    integer, intent(in) :: i, ii, iii, iiii, iiiii
    double precision, dimension(*), intent(in) :: UU
    double precision :: discrete
    double precision, dimension(nw) :: cw
    double precision, dimension(nx) :: cx
    double precision, dimension(ny) :: cy
    double precision, dimension(nz) :: cz
    double precision, dimension(nt) :: ct

    double precision :: w, x, y, z, t, discrete
    integer :: j, jj, jjj, jjjj, jjjjj

    call rderiv1(word, nw, i, hw, cw)
    w = wg(i)
    call rderiv1(xord, nx, ii, hx, cx)
    x = xg(ii)
    call rderiv1(yord, ny, iii, hy, cy)
    y = yg(iii)
    call rderiv1(zord, nz, iiii, hz, cz)
    z = zg(iiii)
    call rderiv1(tord, nt, iiiii, ht, ct)
    t = tg(iiiii)
    discrete = 0.0

    ! five cases of one partial w, x, y, z, t to any order
    if(word.ne.0.and.xord.eq.0.and.yord.eq.0.and.zord.eq.0.and.tord.eq.0) then
      ! Uw  Uww  Uwww  Uwwww
      do j=1,nw
        if(j.eq.1 .or. j.eq.nw) then
          discrete = discrete + cw(j)*uana(wg(j),x,y,z,t)
        else
          discrete = discrete + cw(j)*UU(s(j,ii,iii,iiii,iiiii))
        end if
      end do ! j
    else if(word.eq.0.and.xord.ne.0.and.yord.eq.0.and.zord.eq.0.and.tord.eq.0) then
      ! Ux  Uxx  Uxxx  Uxxxx
      do jj=1,nx
        if(jj.eq.1 .or. jj.eq.nx) then
          discrete = discrete + cx(jj)*uana(w,xg(jj),y,z,t)
        else
          discrete = discrete + cx(jj)*UU(s(i,jj,iii,iiii,iiiii))
        end if
      end do ! jj
    else if(word.eq.0.and.xord.eq.0.and.yord.ne.0.and.zord.eq.0.and.tord.eq.0) then
      ! Uy  Uyy  Uyyy  Uyyyy
      do jjj=1,ny
        if(jjj.eq.1 .or. jjj.eq.ny) then
          discrete = discrete + cy(jjj)*uana(w,x,yg(jjj),z,t)
        else
          discrete = discrete + cy(jj)*UU(s(i,ii,jjj,iiii,iiiii))
        end if
      end do ! jjj
    else if(word.eq.0.and.xord.eq.0.and.yord.eq.0.and.zord.ne.0.and.tord.eq.0) then
      ! Uz  Uzz  Uzzz  Uzzzz
      do jjjj=1,nz
        if(jjjj.eq.1 .or. jjjj.eq.nz) then
          discrete = discrete + cz(jjjj)*uana(w,x,y,zg(jjjj),t)
        else
          discrete = discrete + cz(jjjj)*UU(s(i,ii,iii,jjjj,iiiii))
        end if
      end do ! jjjj
    else if(word.eq.0.and.xord.eq.0.and.yord.eq.0.and.zord.eq.0.and.tord.ne.0) then
      ! Ut  Utt  Uttt  Utttt
      do jjjjj=1,nt
        if(jjjjj.eq.1 .or. jjjjj.eq.nt) then
          discrete = discrete + ct(jjjjj)*uana(w,x,y,z,tg(jjjjj))
        else
           discrete = discrete + ct(jjjjj)*UU(s(i,ii,iii,iiii,jjjjj))
        end if
      end do ! jjjjj
    end if 
    ! return discrete
  end function eval_derivative

  subroutine check_soln(UU) ! for specific PDE
    ! solve Uwwww(w,x,y,z,t) + Uxxxx(w,x,y,z,t) + Uyyyy(w,x,y,z,t) +
    !       Uzzzz(w,x,y,z,t) + Utttt(w,x,y,z,t) + 2 uww(w,x,y,z,t) +
    !       2 Uxx(w,x,y,z,t) + 2 Uyy(w,x,y,z,t) + 2 Uzz(w,x,y,z,t) +
    !       2 Utt(w,x,y,z,t) + 5 U(w,x,y,z,t) = 0

    use global
    implicit none
    double precision, dimension(*), intent(in) :: UU
    double precision :: val, err, smaxerr
    integer :: i, ii, iii, iiii, iiiii

    smaxerr = 0.0
    do i=2,nw-1
      do ii=2,nx-1
        do iii=2,ny-1
          do iiii=2,nz-1
            do iiiii=2,nt-1
              val = 0.0
              ! Uwwww(w,x,y,z,t)
              val = val + eval_derivative(4,0,0,0,0,i,ii,iii,iiii,iiiii,UU)
              ! + Uxxxx(w,x,y,z,t)
              val = val + eval_derivative(0,4,0,0,0,i,ii,iii,iiii,iiiii,UU)
              ! + Uyyyy(w,x,y,z,t)
              val = val + eval_derivative(0,0,4,0,0,i,ii,iii,iiii,iiiii,UU)
              ! + Uzzzz(w,x,y,z,t)
              val = val + eval_derivative(0,0,0,4,0,i,ii,iii,iiii,iiiii,UU)
              ! + Utttt(w,x,y,z,t)
              val = val + eval_derivative(0,0,0,0,4,i,ii,iii,iiii,iiiii,UU)
              ! + 2 uww(w,x,y,z,t)
              val = val + 2.0*eval_derivative(2,0,0,0,0,i,ii,iii,iiii,iiiii,UU)
              ! + 2 uxx(w,x,y,z,t)
              val = val + 2.0*eval_derivative(0,2,0,0,0,i,ii,iii,iiii,iiiii,UU)
              ! + 2 uyy(w,x,y,z,t)
              val = val + 2.0*eval_derivative(0,0,2,0,0,i,ii,iii,iiii,iiiii,UU)
              ! + 2 uzz(w,x,y,z,t)
              val = val + 2.0*eval_derivative(0,0,0,2,0,i,ii,iii,iiii,iiiii,UU)
              ! + 2 utt(w,x,y,z,t)
              val = val + 2.0*eval_derivative(0,0,0,0,2,i,ii,iii,iiii,iiiii,UU)
              ! + 5 u(w,x,y,z,t)
              val = val + 5.0*uana(wg(i),xg(ii),yg(iii),zg(iiii),tg(iiiii))
              ! -       no f(w,x,y,z,t) should be zero
              err = abs(val)
              if(err>smaxerr) then
                smaxerr = err
              end if
            end do !iiiii
          end do !iiii
        end do !iii
      end do ! ii
    end do ! i
    print *, 'check_soln maxerr=',smaxerr
  end subroutine check_soln

  subroutine write_soln(UU) ! for specific PDE
    use global
    implicit none
    double precision, dimension(*), intent(in) :: UU
    integer :: i, ii, iii, iiii, iiiii

    print *, 'write file pde45h_eq_f90.dat'
    open (unit = 11, file = "pde45h_eq_f90.dat", action = "write")
    write(unit = 11, fmt = "(6I3)")&
          5,nw,nx,ny,nz,nt
    do i=2,nx-1
      do ii=2,nx-1
        do iii=2,ny-1
          do iiii=2,nz-1
            do iiiii=2,nt-1
            write(unit = 11, fmt = "(6F12.6)")&
                  wg(i),xg(ii),yg(iii),zg(iiii),tg(iiiii),&
                  UU(s(i,ii,iii,iiii,iiiii))
            end do !iiiii
          end do !iiii
        end do !iii
      end do ! ii
    end do ! i
    print *, 'file written '
  end subroutine write_soln

program pde45h_eq
  use global
  implicit none
  integer :: i, j, ii, jj, iii, jjj, iiii, jjjj, iiiii, jjjjj
  integer :: k
  double precision :: w, x, y, z, t
  double precision :: val, err, avgerr, maxerr

  print *, 'pde45h_eq.f90 running'
  print *,' solve Uwwww(w,x,y,z,t) + Uxxxx(w,x,y,z,t) + Uyyyy(w,x,y,z,t) +'
  print *,'       Uzzzz(w,x,y,z,t) + Utttt(w,x,y,z,t) + 2 Uww(w,x,y,z,t) +'
  print *,'       2 Uxx(w,x,y,z,t) + 2 Uyy(w,x,y,z,t) + 2 Uzz(w,x,y,z,t) +'
  print *,'       2 Utt(w,x,y,z,t) + 2 U(w,x,y,z,t) = 0 '
  print *,' '
  print *,' boundary conditions computed using U(w,x,y,z,t)'
  print *,' analytic solution is U(w,x,y,z,t) = sin(w+x+y+z+t)'
  print *, ' '
  print *, 'wmin=',wmin,', wmax =',wmax
  print *, 'xmin=',xmin,', xmax =',xmax
  print *, 'ymin=',ymin,', ymax =',ymax
  print *, 'zmin=',zmin,', zmax =',zmax
  print *, 'tmin=',tmin,', tmax =',tmax
  print *, 'nw=',nw,', nx=',nx,', ny=',ny,', nz=',nz,', nt=',nt 

  hw = (wmax-wmin)/(nw-1) ! does not have to be uniform spacing
  print *, 'w grid and analytic solution, at xmin,ymin,zmin,tmin'
  do i=1,nw
    wg(i) = wmin+(i-1)*hw
    Ua(i) = uana(wg(i),xmin,ymin,zmin,tmin) ! not saved
    print *, 'i=',i,', Ua(',wg(i),')=',Ua(i)
  end do
  print *, ' '

  hx = (xmax-xmin)/(nx-1) ! does not have to be uniform spacing
  print *, 'x grid and analytic solution, at ymin,zmin,tmin'
  do ii=1,nx 
    xg(ii) = xmin+(ii-1)*hx
    Ua(ii) = uana(wmin,xg(ii),ymin,zmin,tmin) ! not saved
    print *, 'ii=',ii,', Ua(',xg(ii),')=',Ua(ii)
  end do
  print *, ' '

  hy = (ymax-ymin)/(ny-1) ! does not have to be uniform spacing
  print *, 'y grid and analytic solution, at xmin,zmin,tmin'
  do iii=1,ny 
    yg(iii) = ymin+(iii-1)*hy
    Ua(iii) = uana(wmin,xmin,yg(iii),zmin,tmin)
    print *, 'iii=',iii,', Ua(',yg(iii),')=',Ua(iii)
  end do
  print *, ' '

  hz = (zmax-zmin)/(nz-1) ! does not have to be uniform spacing
  print *, 'z grid and analytic solution, at xmin,ymin,tmin'
  do iiii=1,nz 
    zg(iiii) = zmin+(iiii-1)*hz
    Ua(iiii) = uana(wmin,xmin,ymin,zg(iiii),tmin)
    print *, 'iiii=',iiii,', Ua(',zg(iiii),')=',Ua(iiii)
  end do
  print *, ' '

  ht = (tmax-tmin)/(nt-1) ! does not have to be uniform spacing
  print *, 't grid and analytic solution, at xmin,ymin,zmin'
  do iiiii=1,nt 
    tg(iiiii) = tmin+(iiiii-1)*ht
    Ua(iiiii) = uana(wmin,xmin,ymin,zmin,tg(iiiii))
    print *, 'iiiii=',iiiii,', Ua(',tg(iiiii),')=',Ua(iiiii)
  end do
  print *, ' '

  ! put solution in: Ua vector
  do i=2,nw-1
    w = wg(i)
    do ii=2,nx-1
      x = xg(ii)
      do iii=2,ny-1
        y = yg(iii)
        do iiii=2,nz-1
          z = zg(iiii)
          do iiiii=2,nt-1
            t = tg(iiiii)
            Ua(s(i,ii,iii,iiii,iiiii)) = uana(w,x,y,z,t)
          end do ! iiiii
        end do ! iiii
      end do ! iii
    end do ! ii
  end do ! i

  ! initialize k and f  inside boundary
  do i=2,nw-1
    do j=2,nw-1 
      do ii=2,nx-1
        do jj=2,nx-1 
          do iii=2,ny-1 
            do jjj=2,ny-1 
              do iiii=2,nz-1 
                do jjjj=2,nz-1 
                  do iiiii=2,nt-1 
                    do jjjjj=2,nt-1 
                      A(s(i,ii,iii,iiii,iiiii),s(j,jj,jjj,jjjj,jjjjj)) = 0.0
                    end do ! jjjjj
                  end do ! iiiii
                end do ! jjjj
              end do ! iiii
            end do ! jjj
          end do ! iii
        end do ! jj
      end do ! ii
    end do ! j
  end do ! i

  do i=2,nw-1
    do ii=2,nx-1
      do iii=2,ny-1 
        do iiii=2,nz-1 
          do iiiii=2,nt-1 
            R(s(i,ii,iii,iiii,iiiii)) = 0.0
          end do ! iiiii
        end do ! iiii
      end do ! iii
    end do ! ii
  end do ! i

  ! compute entries in stiffness matrix
  call buildA()

  call simeq(nwxyzt, nwxyzt, A, R, U)

  print *, 'check check_soln against solution'
  call check_soln(Ua)
  print *, 'check computed solution'
  call check_soln(U)
  call write_soln(U)

  avgerr = 0.0
  maxerr = 0.0
  print *, 'U computed discrete, Ua analytic, error'
  do i=2,nw-1 
    do ii=2,nx-1 
      do iii=2,ny-1 
        do iiii=2,nz-1 
          do iiiii=2,nt-1 
            err = abs(U(s(i,ii,iii,iiii,iiiii))-Ua(s(i,ii,iii,iiii,iiiii)))
            if(err>maxerr) then
              maxerr = err
            end if
            avgerr = avgerr + err
            print "('U(',I1,',',I1,',',I1,',',I1,',',I1,')=',f7.4,', Ua=',f7.4,', err=',e15.7)", &
                    i,ii,iii,iiii,iiiii,U(s(i,ii,iii,iiii,iiiii)),&
                    Ua(s(i,ii,iii,iiii,iiiii)), &
                    (U(s(i,ii,iii,iiii,iiiii))-Ua(s(i,ii,iii,iiii,iiiii)))
          end do ! iiiii
        end do ! iiii
      end do ! iii
    end do ! ii
  end do ! i
  print *, '           maxerr=',maxerr,', avgerr=',avgerr/(nwxyzt)
  print *, 'end pde45h_eq.f90 '

end program pde45h_eq