! pde44h_eq.f90  homogeneous Biharmonic equation in four dimensions
! solve uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +
!       utttt(x,y,z,t) + 2 uxx(x,y,z,t) + 2 uyy(x,y,z,t) +
!       2  uzz(x,y,z,t) + 2 utt(x,y,z,t) + 4 u(x,y,z,t) = 0
!
! boundary conditions computed using u(x,y,z,t)
! analytic solution is u(x,y,z,t) = sin(x+y+z+t)
!
! solve for Uijkl numerical approximation U(x_i,y_j,z_k,t_l)
! A * U = R, solve simultaneous equations for U

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

  interface
    function uana(x,y,z,t) result(v) ! analytic solution and boundary
      double precision, intent(in) :: 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

  end interface

end module global

module interf
  interface

    function s(i, ii, iii, iiii) result(v)
      integer, intent(in) :: i, ii, iii, iiii
      integer :: v
    end function s

    subroutine buildA()
    end subroutine buildA
    subroutine printA()
    end subroutine printA

    function eval_derivative(xord, yord, zord, tord, &
                             i, ii, iii, iiii, UU) result (discrete)
      integer, intent(in) :: xord, yord, zord, tord, i, ii, iii, iiii
      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 interf

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

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

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

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

  subroutine buildA() ! dpends on specific PDE
    use global
    implicit none
    interface
      function s(i, ii, iii, iiii) result(v)
        integer, intent(in) :: i, ii, iii, iiii
        integer :: v
      end function s
    end interface

    integer :: i, j, ii, jj, iii, jjj, iiii, jjjj, k
    double precision :: val = 0
    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,nx,1,hx,cxxxx)'
    call rderiv1(4,nx,1,hx,cxxxx)
    do i=1,nx 
      print *, 'cxxxx(',i,')=',cxxxx(i)
    end do
    print *, 'check rderiv1(4,nx,nx,hx,cxxxx)'
    call rderiv1(4,nx,nx,hx,cxxxx)
    do i=1,nx 
      print *, 'cxxxx(',i,')=',cxxxx(i)
    end do
    print *, ' '

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

            ! do each term

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

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

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

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

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

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

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

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

            ! do 4 u(x,y,z,t)
            A(k,s(i,ii,iii,iiii)) = A(k,s(i,ii,iii,iiii)) + 4.0

            ! do f(x,y,z,t) RHS, none

            ! end terms

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

  subroutine printA()
    use global
    implicit none
    interface
      function s(i, ii, iii, iiii) result(v)
        integer, intent(in) :: i, ii, iii, iiii
        integer :: v
      end function s
    end interface
    integer :: i, j, ii, jj, iii, jjj, iiii, jjjj, k
    do i=2,nx-1
      do ii=2,ny-1 
        do iii=2,nz-1 
          do iiii=2,nt-1 
            k = s(i,ii,iii,iiii)
            do j=2,nx-1 
              do jj=2,ny-1 
                do jjj=2,nz-1 
                  do jjjj=2,nt-1 
                    print "('A((',I1,',',I1,','I1,',',I1,'),(',I1,',',I1,','I1,',',I1,'))=', e18.10)", &
                          i, ii, iii, iiii, j, jj, jjj, jjjj, &
                          A(k,s(j,jj,jjj,jjjj))
                  end do ! jjjj
                end do ! jjj
              end do ! jj
            end do ! j
            print *, 'RHS(',k,')=',R(k)
          end do ! iiii
        end do ! iii
      end do ! ii
    end do ! i
  end subroutine printA

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

    integer, parameter ::nn =  max(nx, max(ny, max(nz, nt)))
    double precision, dimension(nn) :: p1, p2, p3
    double precision :: x, y, z, t
    integer :: j, jj, jjj, jjjj

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

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

    ! six cases of two partials xy, xz, xt, yz, yt, zt to any order
    else if(xord.ne.0 .and. yord.ne.0 .and. zord.eq.0 .and. tord.eq.0) then
      ! Ux^n y^m
      do j=1,nx
        p1(j) = 0.0
        do jj=1,ny
          if(j.eq.1 .or. j.eq.nx .or. jj.eq.1 .or. jj.eq.ny) then
            p1(j) = p1(j) + cy(jj)*uana(xg(j),yg(jj),zg(iii),tg(iiii))
          else
             p1(j) = p1(j) + cy(jj)*UU(s(j,jj,iii,iiii))
          end if
        end do ! jj
        discrete = discrete + cx(j)*p1(j)
      end do ! j
    else if(xord.ne.0 .and. yord.eq.0 .and. zord.ne.0 .and. tord.eq.0) then
      ! Ux^n z^m
      do j=1,nx
        p1(j) = 0.0
        do jjj=1,nz
          if(j.eq.1 .or. j.eq.nx .or. jjj.eq.1 .or. jjj.eq.nz) then
            p1(j) = p1(j) + ct(jjj)*uana(xg(j),yg(ii),zg(jjj),tg(iiii))
          else
            p1(j) = p1(j) + cz(jjj)*UU(s(j,ii,jjj,iiii))
          end if
        end do ! jjj
        discrete = discrete + cx(j)*p1(j)
      end do ! j
    else if(xord.ne.0 .and. yord.eq.0 .and. zord.eq.0 .and. tord.ne.0) then
      ! Ux^n t^m
      do j=1,nx
        p1(j) = 0.0
        do jjjj=1,nt
          if(j.eq.1 .or. j.eq.nx .or. jjjj.eq.1 .or. jjjj.eq.nt) then
            p1(j) = p1(j) + ct(jjjj)*uana(xg(j),yg(ii),zg(iii),tg(jjjj))
          else
            p1(j) = p1(j) + ct(jjjj)*UU(s(j,ii,iii,jjjj))
          end if
        end do ! jjjj
        discrete = discrete + cx(j)*p1(j)
      end do ! j
    else if(xord.eq.0 .and. yord.ne.0 .and. zord.ne.0 .and. tord.eq.0) then
      ! Uy^n z^m
      do jj=1,ny
        p1(jj) = 0.0
        do jjj=1,nz
          if(jj.eq.1 .or. jj.eq.ny .or. jjj.eq.1 .or. jjj.eq.nz) then
            p1(jj) = p1(jj) + cz(jjj)*uana(xg(i),yg(jj),zg(jjj),tg(iiii))
          else
            p1(jj) = p1(jj) + cz(jjj)*UU(s(i,jj,jjj,iiii))
          end if
        end do ! jjj
        discrete = discrete + cy(jj)*p1(jj)
      end do ! jj
    else if(xord.eq.0 .and. yord.ne.0 .and. zord.eq.0 .and. tord.ne.0) then
      ! Uy^n t^m
      do jj=1,ny
        p1(jj) = 0.0
        do jjjj=1,nt
          if(jj.eq.1 .or. jj.eq.ny .or. jjjj.eq.1 .or. jjjj.eq.nt) then
             p1(jj) = p1(jj) + ct(jjjj)*uana(xg(i),yg(jj),zg(iii),tg(jjjj))
          else
             p1(jj) = p1(jj) + ct(jjjj)*UU(s(i,jj,iii,jjjj))
          end if
        end do ! jjjj
        discrete = discrete + cy(jj)*p1(jj)
      end do ! jj
    else if(xord.eq.0 .and. yord.eq.0 .and. zord.ne.0 .and. tord.ne.0) then
      ! Uz^n t^m
      do jjj=1,nz
        p1(jjj) = 0.0
        do jjjj=1,nt
          if(jjj.eq.1 .or. jjj.eq.nz .or. jjjj.eq.1 .or. jjjj.eq.nt) then
            p1(jjj) = p1(jjj) + ct(jjjj)*uana(xg(i),yg(ii),zg(jjj),tg(jjjj))
          else
            p1(jjj) = p1(jjj) + ct(jjjj)*UU(s(i,ii,jjj,jjjj))
          end if
        end do ! jjjj
        discrete = discrete + cz(jjj)*p1(jjj)
      end do ! jjj

    ! four cases of three partials xyz, xyt, xzt, yzt to any order
    else if(xord.ne.0 .and. yord.ne.0 .and. zord.ne.0 .and. tord.eq.0) then
      ! Ux^n y^m z^p
      do j=1,nx
        p1(j) = 0.0
        do jj=1,ny
          p2(jj) = 0.0
          do jjj=1,nz
            if(j.eq.1 .or. j.eq.nx .or.      &
               jj.eq.1 .or. jj.eq.ny .or.    &
               jjj.eq.1 .or. jjj.eq.nz) then
               p2(jj) = p2(jj) + cz(jjj)*uana(xg(j),yg(jj),zg(jjj),tg(iiii))
            else
               p2(jj) = p2(jj) + cz(jjj)*UU(s(j,jj,jjj,iiii))
            end if
          end do ! jjj
          p1(j) = p1(j) + cy(jj)*p2(jj)
        end do ! jj
        discrete = discrete + cx(j)*p1(j)
      end do ! j
    else if(xord.ne.0 .and. yord.ne.0 .and. zord.eq.0 .and. tord.ne.0) then
      ! Ux^n y^m t^p
      do j=1,nx
        p1(j) = 0.0
        do jj=1,ny
          p2(jj) = 0.0
          do jjjj=1,nt
            if(j.eq.1 .or. j.eq.nx .or.   &
               jj.eq.1 .or. jj.eq.ny .or. &
               jjjj.eq.1 .or. jjjj.eq.nt) then
              p2(jj) = p2(jj) + ct(jjjj)*uana(xg(j),yg(jj),zg(iii),tg(jjjj))
            else
              p2(jj) = p2(jj) + ct(jjjj)*UU(s(j,jj,iii,jjjj))
            end if
          end do ! jjjj
          p1(j) = p1(j) + cy(jj)*p2(jj)
        end do ! jj
        discrete = discrete +cx(j)*p1(j)
      end do ! j
    else if(xord.ne.0 .and. yord.eq.0 .and. zord.ne.0 .and. tord.ne.0) then
      ! Ux^n z^m t^p
      do j=1,nx
        p1(j) = 0.0
        do jjj=1,nz
          p2(jjj) = 0.0
          do jjjj=1,nt
            if(j.eq.1 .or. j.eq.nx .or.      &
               jjj.eq.1 .or. jjj.eq.nz .or.  &
               jjjj.eq.1 .or. jjjj.eq.nt) then
              p2(jjj) = p2(jjj) + ct(jjjj)*uana(xg(j),yg(ii),zg(jjj),tg(jjjj))
            else
              p2(jjj) = p2(jjj) + ct(jjjj)*UU(s(j,ii,jjj,jjjj))
            end if
          end do ! jjjj
          p1(j) = p1(j) + cz(jjj)*p2(jjj)
        end do ! jjj
        discrete = discrete + cx(j)*p1(j)
      end do ! j
    else if(xord.eq.0 .and. yord.ne.0 .and. zord.ne.0 .and. tord.ne.0) then
      ! Uy^n z^m t^p
      do jj=1,ny
        p1(jj) = 0.0
        do jjj=1,nz
          p2(jjj) = 0.0
          do jjjj=1,nt
            if(jj.eq.1 .or. jj.eq.ny .or.   &
               jjj.eq.1 .or. jjj.eq.nz .or. &
               jjjj.eq.1 .or. jjjj.eq.nt) then
              p2(jjj) = p2(jjj) + ct(jjjj)*uana(xg(i),yg(jj),zg(jjj),tg(jjjj))
            else
              p2(jjj) = p2(jjj) + ct(jjjj)*UU(s(i,jj,jjj,jjjj))
            end if
          end do ! jjjj
          p1(jj) = p1(jj) + cz(jjj)*p2(jjj)
        end do ! jjj
        discrete = discrete + cy(jj)*p1(jj)
      end do ! jj

    ! one case of four partials
    else
      ! Uxyzt
      do j=1,nx
        p1(j) = 0.0
        do jj=1,ny
          p2(jj) = 0.0
          do jjj=1,nz
            p3(jjj) = 0.0
            do jjjj=1,nt
              if(j.eq.1 .or. j.eq.nx .or.     &
                 jj.eq.1 .or. jj.eq.ny .or.   &
                 jjj.eq.1 .or. jjj.eq.nz .or. &
                 jjjj.eq.1 .or. jjjj.eq.nt) then
                p3(jjj) = p3(jjj) + ct(jjjj)*uana(xg(j),yg(jj),zg(jjj),tg(jjjj))
              else
                p3(jjj) =  p3(jjj) + ct(jjjj)*UU(s(j,jj,jjj,jjjj))
              end if
            end do ! jjjj
            p2(jj) = p2(jj) + cz(jjj)*p3(jjj)
          end do ! jjj
          p1(j) = p1(j) + cy(jj)*p2(jj)
        end do ! jj
        discrete = discrete + cx(j)*p1(j)
      end do ! j
    end if
    ! return discrete
  end function eval_derivative

  subroutine check_soln(UU) ! for specific PDE
    ! solve uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +
    !       utttt(x,y,z,t) + 2 uxx(x,y,z,t) + 2 uyy(x,y,z,t) +
    !       2 uzz(x,y,z,t) + 2 utt(x,y,z,t) + 2 u(x,y,z,t) = 0

    use global
    implicit none
    interface
      function s(i, ii, iii, iiii) result(v)
        integer, intent(in) :: i, ii, iii, iiii
        integer :: v
      end function s
      function eval_derivative(xord, yord, zord, tord, &
                               i, ii, iii, iiii, UU) result (discrete)
        integer, intent(in) :: xord, yord, zord, tord, i, ii, iii, iiii
        double precision, dimension(*), intent(in) :: UU
        double precision :: discrete
      end function eval_derivative
    end interface
    double precision, dimension(*), intent(in) :: UU
    double precision :: val, err, smaxerr
    integer :: i, ii, iii, iiii

    smaxerr = 0.0
    do i=2,nx-1
      do ii=2,ny-1
        do iii=2,nz-1
          do iiii=2,nt-1
            val = 0.0
            ! uxxxx(x,y,z,t)
            val = val + eval_derivative(4, 0, 0, 0, i, ii, iii, iiii, UU)
            ! + uyyyy(x,y,z,t)
            val = val + eval_derivative(0, 4, 0, 0, i, ii, iii, iiii, UU)
            ! + uzzzz(x,y,z,t)
            val = val + eval_derivative(0, 0, 4, 0, i, ii, iii, iiii, UU)
            ! + utttt(x,y,z,t)
            val = val + eval_derivative(0, 0, 0, 4, i, ii, iii, iiii, UU)
            ! + 2  uxx(x,y,z,t)
            val = val + 2.0*eval_derivative(2, 0, 0, 0, i, ii, iii, iiii, UU)
            ! + 2 uyy(x,y,z,t)
            val = val + 2.0*eval_derivative(0, 2, 0, 0, i, ii, iii, iiii, UU)
            ! + 2  uzz(x,y,z,t)
            val = val + 2.0*eval_derivative(0, 0, 2, 0, i, ii, iii, iiii, UU)
            ! + 2  utt(x,y,z,t)
            val = val + 2.0*eval_derivative(0, 0, 0, 2, i, ii, iii, iiii, UU)
            ! + 4 u(x,y,z,t)
            val = val + 4.0*uana(xg(i),yg(ii),zg(iii),tg(iiii))
            ! -       no f(x,y,z,t) should be zero
            err = abs(val)
            if(err>smaxerr) then
              smaxerr = err
            end if
          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
    interface
      function s(i, ii, iii, iiii) result(v)
        integer, intent(in) :: i, ii, iii, iiii
        integer :: v
      end function s
    end interface
    double precision, dimension(*), intent(in) :: UU
    integer :: i, ii, iii, iiii

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

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

  print *, 'pde44h_eq.f90 running'
  print *,' solve uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +'
  print *,'       utttt(x,y,z,t) + 2 uxx(x,y,z,t) + 2 uyy(x,y,z,t) +'
  print *,'       2 uzz(x,y,z,t) + 2 utt(x,y,z,t) + 4 u(x,y,z,t) = 0'
  print *,' '
  print *,' boundary conditions computed using u(x,y,z,t)'
  print *,' analytic solution is u(x,y,z,t) = sin(x+y+z+t)'
  print *, ' '
  print *, 'xmin=',xmin,', xmax =',xmax
  print *, 'ymin=',ymin,', ymax =',ymax
  print *, 'zmin=',zmin,', zmax =',zmax
  print *, 'tmin=',tmin,', tmax =',tmax
  print *, 'nx=',nx,', ny =',ny,', nz =',nz,', nt=',nt 

  hx = (xmax-xmin)/(nx-1) ! does not have to be uniform spacing
  print *, 'x grid and analytic solution, at ymin,zmin,tmin'
  do i=1,nx 
    xg(i) = xmin+(i-1)*hx
    Ua(i) = uana(xg(i),ymin,zmin,tmin) ! not saved
    print *, 'i=',i,', Ua(',xg(i),')=',Ua(i)
  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 ii=1,ny 
    yg(ii) = ymin+(ii-1)*hy
    Ua(ii) = uana(xmin,yg(ii),zmin,tmin)
    print *, 'ii=',ii,', Ua(',yg(ii),')=',Ua(ii)
  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 iii=1,nz 
    zg(iii) = zmin+(iii-1)*hz
    Ua(iii) = uana(xmin,ymin,zg(iii),tmin)
    print *, 'iii=',iii,', Ua(',zg(iii),')=',Ua(iii)
  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 iiii=1,nt 
    tg(iiii) = tmin+(iiii-1)*ht
    Ua(iiii) = uana(xmin,ymin,zmin,tg(iiii))
    print *, 'iiii=',iiii,', Ua(',tg(iiii),')=',Ua(iiii)
  end do
  print *, ' '

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

  ! initialize k and f  inside boundary
  do i=2,nx-1
    do j=2,nx-1 
      do ii=2,ny-1 
        do jj=2,ny-1 
          do iii=2,nz-1 
            do jjj=2,nz-1 
              do iiii=2,nt-1 
                do jjjj=2,nt-1 
                  A(s(i,ii,iii,iiii),s(j,jj,jjj,jjjj)) = 0.0
                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,nx-1
    do ii=2,ny-1 
      do iii=2,nz-1 
        do iiii=2,nt-1 
          R(s(i,ii,iii,iiii)) = 0.0
        end do ! iiii
      end do ! iii
    end do ! ii
  end do ! i

  ! compute entries in stiffness matrix
  call buildA()

  print *, 'A and R computed'
  ! call printA() ! optional
  print *, ' '

  call simeq(nxyzt, nxyzt, 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 Galerkin, Ua analytic, error'
  do i=2,nx-1 
    do ii=2,ny-1 
      do iii=2,nz-1 
        do iiii=2,nt-1 
          err = abs(U(s(i,ii,iii,iiii))-Ua(s(i,ii,iii,iiii)))
          if(err>maxerr) then
            maxerr = err
          end if
          avgerr = avgerr + err
          print "('U(',I1,',',I1,','I1,',',I1')=',f7.4,', Ua=',f7.4,', err=',e15.7)", &
              i, ii, iii, iiii, U(s(i,ii,iii,iiii)), Ua(s(i,ii,iii,iiii)), &
              (U(s(i,ii,iii,iiii))-Ua(s(i,ii,iii,iiii)))
        end do ! iiii
      end do ! iii
    end do ! ii
  end do ! i
  print *, '           maxerr=',maxerr,', avgerr=',avgerr/(nxyzt)
  print *, 'end pde44h_eq.f90 '

end program pde44h_eq