! pde45h_eq.f90  homogeneous Biharmonic equation in five dimensions
! solve Uxxxx(x,y,z,t,u) + Uyyyy(x,y,z,t,u) + Uzzzz(x,y,z,t,u) + 
!       Utttt(x,y,z,t,u) + Uuuuu(x,y,z,t,u) +
!       2 Uxx(x,y,z,t,u) + 2 Uyy(x,y,z,t,u) + 2 Uzz(x,y,z,t,u) +
!       2 Utt(x,y,z,t,u) + 2 Uuu(x,y,z,t,u) + 5 U(x,y,z,t,u) = 0
!
! boundary conditions computed using U(x,y,z,t,u)
! analytic solution is U(x,y,z,t,u) = sin(x+y+z+t+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 :: nu=6
  integer, parameter :: nxyztu=(nx-2)*(ny-2)*(nz-2)*(nt-2)*(nu-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 :: umin = 0.0
  double precision :: umax = 0.5
  double precision :: hx, hy, hz, ht, hu
  double precision, dimension(nxyztu,nxyztu) :: A ! stiffness matrix
     ! A(s(i,ii,i3,i4,i5),s(j,jj,j3,j4,j5))
  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(nu) :: ug ! u
  double precision, dimension(nxyztu) :: R ! RHS
  double precision, dimension(nxyztu) :: Uk ! unknown solution
  double precision, dimension(nxyztu) :: Ua ! known solution for checking

  interface
    function uana(x,y,z,t,u) result(v) ! analytic solution and boundary
      double precision, intent(in) :: x,y,z,t,u
      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 rderiv(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 rderiv
    subroutine rderiv1(order, npoints, point, hx, cx) ! point-1
      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 formain
  interface
    subroutine buildA()
    end subroutine buildA
    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 formain

module uses
  interface
    function s(i, ii, i3, i4, i5) result(v)
      integer, intent(in) :: i, ii, i3, i4, i5
      integer :: v
    end function s
  end interface
end module uses

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

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

  function s(i, ii, i3, i4, i5) result(v)
    use global
    implicit none
    integer, intent(in) :: i, ii, i3, i4, i5
    integer :: v

    v = (i-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2) + &
        (ii-2)*(ny-2)*(nz-2)*(nt-2) + (i3-2)*(nz-2)*(nt-2) + &
        (i4-2)*(nt-2) + i5-1
  end function s

  subroutine buildA() ! dpends on specific PDE
    use global
    use uses
    implicit none
    integer :: i, j, ii, jj, i3, j3, i4, j4, i5, j5, 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
    double precision, dimension(nu) :: cuu, cuuuu

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


    print *, 'check nuderiv(4,nx,1,hx,cxxxx)'
    call nuderiv(4,nx,1,xg,cxxxx)
    do i=1,nx 
      print *, 'cxxxx(',i,')=',cxxxx(i)
    end do
    print *, 'check rderiv(4,nx,0,hx,cxxxx)'
    call rderiv(4,nx,0,hx,cxxxx)
    do i=1,nx 
      print *, 'cxxxx(',i,')=',cxxxx(i)
    end do
    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,nu,1,hu,cuuuu)'
    call rderiv1(4,nu,1,hu,cuuuu)
    do i=1,nu 
      print *, 'cuuuu(',i,')=',cuuuu(i)
    end do
    print *, 'check rderiv1(4,nu,nu,hu,cuuuu)'
    call rderiv1(4,nu,nu,hu,cuuuu)
    do i=1,nu
      print *, 'cuuuu(',i,')=',cuuuu(i)
    end do
    print *, ' '

    val = 0.0
    do i=2,nx-1
      do ii=2,ny-1
        do i3=2,nz-1
          do i4=2,nt-1
            do i5=2,nu-1
              k = s(i,ii,i3,i4,i5) ! 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(i3),tg(i4),ug(i5))
                else
                  A(k,s(j,ii,i3,i4,i5))=A(k,s(j,ii,i3,i4,i5)) + 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(i3),tg(i4),ug(i5))
                else
                  A(k,s(i,jj,i3,i4,i5))=A(k,s(i,jj,i3,i4,i5)) + val
                end if
              end do ! jj

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

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

              ! do d^4U/du^4
              call rderiv1(4,nu,i5,hu,cuuuu)
              do j5=1,nu
                val = cuuuu(j5)
                if(j5.eq.1 .or. j5.eq.nu) then
                  R(k) = R(k)-val*uana(xg(i),yg(ii),zg(i3),tg(i4),ug(j5))
                else
                  A(k,s(i,ii,i3,i4,j5))=A(k,s(i,ii,i3,i4,j5)) + val
                end if
              end do ! j5

              ! 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(i3),tg(i4),ug(i5))
                else
                  A(k,s(j,ii,i3,i4,i5))=A(k,s(j,ii,i3,i4,i5)) + 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(i3),tg(i4),ug(i5))
                else
                  A(k,s(i,jj,i3,i4,i5))=A(k,s(i,jj,i3,i4,i5)) + val
                end if
              end do ! jj

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

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

              ! do 2 d^2U/du^2
              call rderiv1(2,nu,i5,hu,cuu)
              do j5=1,nu
                val = 2.0*cuu(j5)
                if(j5.eq.1 .or. j5.eq.nu) then
                  R(k) = R(k)-val*uana(xg(i),yg(ii),zg(i3),tg(i4),ug(j5))
                else
                  A(k,s(i,ii,i3,i4,j5))=A(k,s(i,ii,i3,i4,j5)) + val
                end if
              end do ! j5

              ! do 5 u(x,y,z,t,u)
              A(k,s(i,ii,i3,i4,i5))=A(k,s(i,ii,i3,i4,i5)) + 5.0

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

              ! end terms

            end do ! i5
          end do ! i4
        end do ! i3
      end do ! ii 
    end do ! i
  end subroutine buildA

  function eval_derivative(xord, yord, zord, tord, uord, &
                           i, ii, i3, i4, i5, UU) result (discrete)
    use global
    use uses
    implicit none
    integer, intent(in) :: xord, yord, zord, tord, uord
    integer, intent(in) :: i, ii, i3, i4, i5
    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
    double precision, dimension(nu) :: cu

    double precision :: x, y, z, t, u
    integer :: j, jj, j3, j4, j5

    call rderiv1(xord, nx, i, hx, cx)
    x = xg(i)
    call rderiv1(yord, ny, ii, hy, cy)
    y = yg(ii)
    call rderiv1(zord, nz, i3, hz, cz)
    z = zg(i3)
    call rderiv1(tord, nt, i4, ht, ct)
    t = tg(i4)
    call rderiv1(uord, nu, i5, hu, cu)
    u = ug(i5)
    discrete = 0.0

    ! five cases of one partial  x, y, z, t, u  to any order
    if(xord.ne.0 .and. yord.eq.0 .and. zord.eq.0 .and. tord.eq.0 &
       .and. uord.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,u)
        else
          discrete = discrete + cx(j)*UU(s(j,ii,i3,i4,i5))
        end if
      end do ! j

    else if(xord.eq.0 .and. yord.ne.0 .and. zord.eq.0 .and. tord.eq.0 &
            .and. uord.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,u)
        else
          discrete = discrete + cy(jj)*UU(s(i,jj,i3,i4,i5))
        end if
      end do ! jj

    else if(xord.eq.0 .and. yord.eq.0 .and. zord.ne.0 .and. tord.eq.0 &
            .and. uord.eq.0) then
      ! Uz  Uzz  Uzzz  Uzzzz
      do j3=1,nz
        if(j3.eq.1 .or. j3.eq.nz) then
          discrete = discrete + cz(j3)*uana(x,y,zg(j3),t,u)
        else
          discrete = discrete + cz(j3)*UU(s(i,ii,j3,i4,i5))
        end if
      end do ! j3

    else if(xord.eq.0 .and. yord.eq.0 .and. zord.eq.0 .and. tord.ne.0 &
            .and. uord.eq.0) then
      ! Ut  Utt  Uttt  Utttt
      do j4=1,nt
        if(j4.eq.1 .or. j4.eq.nt) then
          discrete = discrete + ct(j4)*uana(x,y,z,ug(j4),u)
        else
           discrete = discrete + ct(j4)*UU(s(i,ii,i3,j4,i5))
        end if
      end do ! j4

    else if(xord.eq.0 .and. yord.eq.0 .and. zord.eq.0 .and. tord.eq.0 &
            .and. uord.ne.0) then
      ! Uu  Uuu  Uuuu  Uuuuu
      do j5=1,nu
        if(j5.eq.1 .or. j5.eq.nu) then
          discrete = discrete + cu(j5)*uana(x,y,z,t,ug(j5))
        else
           discrete = discrete + cu(j5)*UU(s(i,ii,i3,i4,j5))
        end if
      end do ! j5

    end if 
    ! return discrete
  end function eval_derivative

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

    use global
    implicit none
    interface
      function eval_derivative(xord, yord, zord, tord, uord, &
                               i, ii, i3, i4, i5, UU) result (discrete)
        integer, intent(in) :: xord, yord, zord, tord, uord
        integer, intent(in) :: i, ii, i3, i4, i5
        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, i3, i4, i5

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

  subroutine write_soln(UU) ! for specific PDE
    use global
    use uses
    implicit none
    double precision, dimension(*), intent(in) :: UU
    integer :: i, ii, i3, i4, i5

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

program pde45h_eq
  use global
  use formain
  use uses
  implicit none
  integer :: i, j, ii, jj, i3, j3, i4, j4, i5, j5
  integer :: k
  double precision :: x, y, z, t, u
  double precision :: val, err, avgerr, maxerr

  print *, 'pde45h_eq.f90 running'
  print *,' solve Uxxxx(x,y,z,t,u) + Uyyyy(x,y,z,t,u) + Uzzzz(x,y,z,t,u) +'
  print *,'       Utttt(x,y,z,t,u) + Uuuuu(x,y,z,t,u) + '
  print *,'       2 Uxx(x,y,z,t,u) + 2 Uyy(x,y,z,t,u) + 2 Uzz(x,y,z,t,u) +'
  print *,'       2 Utt(x,y,z,t,u) + 2 Uuu(x,y,z,t,u) + 5 U(x,y,z,t,u) = 0 '
  print *,' '
  print *,' boundary conditions computed using U(x,y,z,t,u)'
  print *,' analytic solution is U(x,y,z,t,u) = sin(x+y+z+t+u)'
  print *, ' '
  print *, 'xmin=',xmin,', xmax =',xmax
  print *, 'ymin=',ymin,', ymax =',ymax
  print *, 'zmin=',zmin,', zmax =',zmax
  print *, 'tmin=',tmin,', tmax =',tmax
  print *, 'umin=',umin,', umax =',umax
  print *, 'nx=',nx,', ny=',ny,', nz=',nz,', nt=',nt,', nu=',nu 

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

  hu = (umax-umin)/(nu-1) ! does not have to be uniform spacing
  print *, 'u grid and analytic solution, at xmin,ymin,zmin,tmin'
  do i5=1,nu
    ug(i5) = umin+(i5-1)*hu
    Ua(i5) = uana(xmin,ymin,zmin,tmin,ug(i5)) ! not saved
    print *, 'i5=',i5,', Ua(',ug(i5),')=',Ua(i5)
  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 i3=2,nz-1
        z = zg(i3)
        do i4=2,nt-1
          t = tg(i4)
          do i5=2,nu-1
            u = ug(i5)
            Ua(s(i,ii,i3,i4,i5)) = uana(x,y,z,t,u)
          end do ! i5
        end do ! i4
      end do ! i3
    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 i3=2,nz-1 
            do j3=2,nz-1 
              do i4=2,nt-1 
                do j4=2,nt-1 
                  do i5=2,nu-1 
                    do j5=2,nu-1 
                      A(s(i,ii,i3,i4,i5),s(j,jj,j3,j4,j5)) = 0.0
                    end do ! j5
                  end do ! i5
                end do ! j4
              end do ! i4
            end do ! j3
          end do ! i3
        end do ! jj
      end do ! ii
    end do ! j
  end do ! i

  do i=2,nx-1
    do ii=2,ny-1
      do i3=2,nz-1 
        do i4=2,nt-1 
          do i5=2,nu-1 
            R(s(i,ii,i3,i4,i5)) = 0.0
          end do ! i5
        end do ! i4
      end do ! i3
    end do ! ii
  end do ! i

  print *, 'compute entries in stiffness matrix'
  call buildA()

  do i=nx-1,nx-1
    do ii=ny-1,ny-1
      do i3=nz-1,nz-1 
        do i4=nt-1,nt-1 
          do i5=nu-1,nu-1 
            k = s(i,ii,i3,i4,i5)
            do j5=2,nu-1
              print *, 'k=',k,' j5=',j5,' A=',A(k,j5)
            end do ! j5
          end do ! i5
        end do ! i4
      end do ! i3
    end do ! ii
  end do ! i

  print *, 'solve ',nxyztu,' equations in ',nxyztu,' unknowns'
  call simeq(nxyztu, nxyztu, A, R, Uk)

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

  avgerr = 0.0
  maxerr = 0.0
  print *, 'U computed discrete, Ua analytic, error'
  do i=2,nx-1 
    do ii=2,ny-1 
      do i3=2,nz-1 
        do i4=2,nt-1 
          do i5=2,nu-1 
            err = abs(Uk(s(i,ii,i3,i4,i5))-Ua(s(i,ii,i3,i4,i5)))
            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,i3,i4,i5,Uk(s(i,ii,i3,i4,i5)),&
                    Ua(s(i,ii,i3,i4,i5)), &
                    (Uk(s(i,ii,i3,i4,i5))-Ua(s(i,ii,i3,i4,i5)))
          end do ! i5
        end do ! i4
      end do ! i3
    end do ! ii
  end do ! i
  print *, '           maxerr=',maxerr,', avgerr=',avgerr/(nxyztu)
  print *, 'end pde45h_eq.f90 '

end program pde45h_eq