// Pde44_eq.scala  PDE 4th order, 4 dimensions
// solve uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +
//     4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +
//     7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +
//     10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) =
//     F(x,y,z,t)
//     F(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +
//                  180 z^2*t + 200 y^2z*t + 44 t^3 +
//                  110 y^2z^2t + 66 xy + 12t^4 +
//                  24 z^4 + 36 y^4 + 48 x^4 +
//                  60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t 
//
// boundary conditions computed using u(x,y,z,t)
// analytic solution is u(x,y,z,t) = t^4 + 2 z^4 + 3 y^4 + 4 x^4 +
//                                   5 y^2 z^2 t^2 + 6 x y t + 
//                                   7 x z + 8 t + 9
//
// replace continuous derivatives with discrete derivatives
// solve for Uijkl numerical approximation U(x_i,y_j,z_k,t_l)
// A * U = F, solve simultaneous equations for U
//
// subscripts in code use x(i), y(ii), z(iii), t(iiii)
// derivatives are computed with b,hx,a included as in
// cx(j) for first derivative of x at point j
// cxxxx(j) for fourth derivative of x at point j
// cyy(jj)*czz(jjj) for d^4u/dy^2dz^2 at jj, jjj
//
// The subscript versions of the derivatives are substituted into the
// equation to be solved. The resulting equation is analytically solved
// for u(i,j,k,l) = other u terms with coefficients and + f(xi,yj,zk,tl)
// The boundaries xmin..xmax, ymin..ymax, zmin..zmax, tmin..tmax are known
// The number of grid points in each dimension is known nx, ny, nz, nt
// hx=(xmax-xmin)/(nx-1), hy=(ymax-ymin)/(ny-1),
// hz=(zmax-zmin)/(nz-1), ht=(tmax-tmin)/(tx-1), 
// thus xi=xmin+i*hx, yj=ymin+j*hy, zk=zmin+k*hz, tl=tmin+l*ht
// then the linear system of equations is built,
// nxyzt=(nx-2)*(ny-2)*(nz-2)*(nt-2) equations because the boundary value
//                                   is known on each end.
// The matrix A is nxyzt rows by nxyzt columns (nxyzt+1) columns including F
// The forcing function vector F is nxyzt elements adjunct to A
// The solution vector U is nxyzt elements.
// The building of the A matix requires 4 nested loops for rows
// i  in 1..nx-1, j  in 1..ny-1,  k  in 1..nz-1, l  in 1..nt-1 for rows of A
// then each term in the differential equation fills in that row
// and cs = (nx-2)*(ny-2)*(nz-2)*(nt-2) the RHS, fills the last column of A
// subscripting functions are used to compute linear subscripts of A, U
// row=(i-1)*(ny-2)*(nz-2)*(nt-2)+(ii-1)*(nz-2)*(nt-2)+(iii-1)*(nt-2)+(iiii-1)
// row*(nxyzt+1)+col for A
//
// care must be taken to move the negative of boundary elements to
// the right hand size of the equation to be solved. These are subtracted
// from the computed value of the forcing function, cs, for that row.
// tests are j==0 || j==nx-1   ...  jjjj==0 || jjjj==nt-1 for boundaries
//

import scala.math._

object Pde44
{
  var S = Simeq
  var D = Deriv

  var nx = 5
  var ny = 5
  var nz = 5
  var nt = 5
  var nxyzt = (nx-2)*(ny-2)*(nz-2)*(nt-2)
  var xmin:Double = 0.0
  var xmax:Double = 1.0
  var ymin:Double = 0.0
  var ymax:Double = 1.0
  var zmin:Double = 0.0
  var zmax:Double = 1.0
  var tmin:Double = 0.0
  var tmax:Double = 1.0
  var xg = Array.ofDim[Double](nx) // grid
  var yg = Array.ofDim[Double](ny)
  var zg = Array.ofDim[Double](nz)
  var tg = Array.ofDim[Double](nt)
  var hx:Double = 0.0
  var hy:Double = 0.0
  var hz:Double = 0.0
  var ht:Double = 0.0
  var A = Array.ofDim[Double](nxyzt,nxyzt+1)

  def main(args: Array[String])
  {
    var U = Array.ofDim[Double](nxyzt)
    var Ua = Array.ofDim[Double](nxyzt)
    var F = Array.ofDim[Double](nxyzt)
    // coeficients of terms in first equation 
    var coefdxxxx:Double = 1.0 // computed below if not constant 
    var coefdyyyy:Double = 2.0
    var coefdzzzz:Double = 3.0
    var coefdtttt:Double = 4.0
    var coefdxyt:Double  = 5.0
    var coefdyyzz:Double = 6.0
    var coefdztt:Double  = 7.0
    var coefdxxx:Double  = 8.0
    var coefdyyt:Double  = 9.0
    var coefdzt:Double   = 10.0
    var coefdt:Double    = 11.0
    var coefu:Double     = 12.0
    // derivative point arrays, dynamic in this version 
    var cx = Array.ofDim[Double](nx) // first derivative
    var cy = Array.ofDim[Double](ny)
    var cz = Array.ofDim[Double](nz)
    var ct = Array.ofDim[Double](nt)
    var cxx = Array.ofDim[Double](nx) // second derivative
    var cyy = Array.ofDim[Double](ny)
    var czz = Array.ofDim[Double](nz)
    var ctt = Array.ofDim[Double](nt)
    var cxxx = Array.ofDim[Double](nx) // third derivative
    var cyyy = Array.ofDim[Double](ny)
    var czzz = Array.ofDim[Double](nz)
    var cttt = Array.ofDim[Double](nt)
    var cxxxx = Array.ofDim[Double](nx) // fourth derivative
    var cyyyy= Array.ofDim[Double](ny)
    var czzzz = Array.ofDim[Double](nz)
    var ctttt = Array.ofDim[Double](nt)
    var err:Double = 0.0
    var avgerr:Double = 0.0
    var maxerr:Double = 0.0
    var valu:Double = 0.0
    var x:Double = 0.0
    var y:Double = 0.0
    var z:Double = 0.0
    var t:Double = 0.0
    var k:Int = 0
    var cs:Int = 0
    var Tmsec = System.currentTimeMillis

    println("Pde44_eq.c running ")
    println("Given uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) + ")
    println("    4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) + ")
    println("    7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) + ")
    println("    10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) = ")
    println("  F(x,y,z,t) ")
    println("  F(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +   ")
    println("               180 z^2*t + 200 y^2z*t + 44 t^3 +    ")
    println("               110 y^2z^2t + 66 xy + 12t^4 +        ")
    println("               24 z^4 + 36 y^4 + 48 x^4 +           ")
    println("               60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t ") 
    println("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries ")
    println("Analytic solution u(x,y,z,t) =  t^4 + 2 z^4 + 3 y^4 + ")
    println("     4 x^4 + 5 y^2 z^2 t^2 + 6 x y t + 7 x z + 8 t + 9 ")
    println("")

    println("xmin="+xmin+", xmax="+xmax+", ymin="+ymin+", ymax="+ymax)
    println("zmin="+zmin+", zmax="+zmax+", tmin="+tmin+", tmax="+tmax)
    println("nx="+nx+", ny="+ny+", nz="+nz+", nt="+nt)

    println("x grid and analytic solution at ymin,zmin,tmin ")
    hx = (xmax-xmin)/((nx-1).toDouble)
    for(i <- 0 until nx )
    {
      xg(i) = xmin+(i.toDouble)*hx
      Ua(i) = uana(xg(i),ymin,zmin,tmin) // temp 
      println("i="+i+", Ua("+xg(i)+","+ymin+","+zmin+","+tmin+")="+Ua(i))
    } // i   
    println("")
    println("y grid and analytic solution at xmin,zmin,tmin ")
    hy = (ymax-ymin)/((ny-1).toDouble)
    for(ii <- 0 until ny )
    {
      yg(ii) = ymin+(ii.toDouble)*hy
      Ua(ii) = uana(xmin, yg(ii),zmin,tmin) // temp 
      println("ii="+ii+", Ua("+xmin+","+yg(ii)+","+zmin+","+tmin+")="+Ua(ii))
    } // ii 
    println("")
    println("z grid and analytic solution at xmin,ymin,tmin ")
    hz = (zmax-zmin)/((nz-1).toDouble)
    for(iii <- 0 until nz )
    {
      zg(iii) = zmin+(iii.toDouble)*hz
      Ua(iii) = uana(xmin, ymin, zg(iii),tmin) // temp 
      println("iii="+iii+", Ua("+xmin+","+ymin+","+zg(iii)+","+tmin+")="+Ua(iii))
    } // iii 
    println("")
    println("t grid and analytic solution at xmin,ymin,zmin ")
    ht = (tmax-tmin)/((nt-1).toDouble)
    for(iiii <- 0 until nt )
    {
      tg(iiii) = tmin+(iiii.toDouble)*ht
      Ua(iiii) = uana(xmin, ymin, zmin, tg(iiii)) // temp 
      println("iiii="+iiii+", Ua("+xmin+","+ymin+","+zmin+","+tg(iiii)+")="+Ua(iiii))
    } // iiii 
    println("")

    println("check rderiv(4,nx,0,hx,cxxxx ")
    cxxxx = D.rderiv(4,nx,0,hx)
    for(i <- 0 until nx )
    {
      println("cxxxx("+i+")="+cxxxx(i))
    }
    println("check rderiv(4,nx,0,hx,cxxxx ")
    cxxxx = D.rderiv(4,nx,nx-1,hx)
    for(i <- 0 until nx )
    {
      println("cxxxx("+i+")="+cxxxx(i))
    }
    println("")

    println("put solution in Ua vector ")
    // put solution in Ua vector 
    for(i <- 1 until nx-1 )
    {
      x = xg(i)
      for(ii <- 1 until ny-1 )
      {
        y = yg(ii)
        for(iii <- 1 until nz-1 )
        {
          z = zg(iii)
          for(iiii <- 1 until nt-1 )
          {
            t = tg(iiii)
            Ua(s(i,ii,iii,iiii)) = uana(x,y,z,t)
          } // iiii 
        } // iii 
      } // ii 
    } // i 
    println("")

    println("initialize A ")
    // initialize A 
    for(i <- 1 until nx-1 )
    {
      for(ii <- 1 until ny-1 )
      {
        for(iii <- 1 until nz-1 )
        {
          for(iiii <- 1 until nt-1 )
          {
            for(j <- 1 until nx-1 )
            {
              for(jj <- 1 until ny-1 )
              {
                for(jjj <- 1 until nz-1 )
                {
                  for(jjjj <- 1 until nt-1 )
                  {
                    A(s(i,ii,iii,iiii))(s(j,jj,jjj,jjjj)) = 0.0
	          } // jjjj 
	        } // jjj 
	      } // jj 
            } // j 
          } // iiii 
        } // iii 
      } // ii 
    } // i 

    // compute entries in A matrix, inside boundary 
    println("compute A matrix ")
    cs = (nx-2)*(ny-2)*(nz-2)*(nt-2) // constant RHS column 
    valu = 0.0
    for(i <- 1 until nx-1 )
    {
      for(ii <- 1 until ny-1 )
      {
        for(iii <- 1 until nz-1 )
        {
          for(iiii <- 1 until nt-1 )
          {
            k = s(i,ii,iii,iiii) // row index 

            // for each term in first equation 

	    // for d^4u/dx^4 
            coefdxxxx = 1.0
            cxxxx = D.rderiv(4,nx,i,hx)
            for(j <- 0 until nx )
            {
              valu = coefdxxxx*cxxxx(j)
              if(j==0 || j==nx-1)
              {
	        A(k)(cs) = A(k)(cs) - valu*uana(xg(j),yg(ii),zg(iii),tg(iiii))
              }
              else
              {
                A(k)(s(j,ii,iii,iiii)) = valu + A(k)(s(j,ii,iii,iiii))
              }
	    }

	    // for d^4u/dy^4 
            coefdyyyy = 2.0
            cyyyy = D.rderiv(4,ny,ii,hy)
            for(jj <- 0 until ny )
            {
              valu = coefdyyyy*cyyyy(jj)
              if(jj==0 || jj==ny-1)
              {
	        A(k)(cs) = A(k)(cs) - valu*uana(xg(i),yg(jj),zg(iii),tg(iiii))
              }
              else
              { 
                A(k)(s(i,jj,iii,iiii)) = valu +A(k)(s(i,jj,iii,iiii))
              }
	    }

	    // for d^4u/dz^4 
            coefdzzzz = 3.0
            czzzz = D.rderiv(4,nz,iii,hz)
            for(jjj <- 0 until nz )
            {
              valu = coefdzzzz*czzzz(jjj)
              if(jjj==0 || jjj==nz-1)
              {
	        A(k)(cs) = A(k)(cs) - valu*uana(xg(i),yg(ii),zg(jjj),tg(iiii))
              }
              else
              {
                A(k)(s(i,ii,jjj,iiii)) = valu + A(k)(s(i,ii,jjj,iiii))
              }
	    }

	    // for d^4u/dt^4 
            coefdtttt = 4.0
            ctttt = D.rderiv(4,nt,iiii,ht)
            for(jjjj <- 0 until nt )
            {
              valu = coefdtttt*ctttt(jjjj)
              if(jjjj==0 || jjjj==nt-1)
              {
	        A(k)(cs) = A(k)(cs) - valu*uana(xg(i),yg(ii),zg(iii),tg(jjjj))
              }
              else
              {
                A(k)(s(i,ii,iii,jjjj)) = valu +A(k)(s(i,ii,iii,jjjj))
              }
	    }

	    // for d^3u/dxdydt 
            coefdxyt  = 5.0
            cx = D.rderiv(1,nx,i,hx)
            cy = D.rderiv(1,ny,ii,hy)
            ct = D.rderiv(1,nt,iiii,ht)
            for(j <- 0 until nx )
	    {
              for(jj <- 0 until ny )
	      {
                for(jjjj <- 0 until nt )
                {
                  valu = coefdxyt*cx(j)*cy(jj)*ct(jjjj)
                  if(j==0 || j==nx-1 || jj==0 || jj==ny-1 ||
                     jjjj==0 || jjjj==nt-1)
                  {
	            A(k)(cs) = A(k)(cs) - valu*uana(xg(j),yg(jj),zg(iii),tg(jjjj))
                  }
                  else
                  {
                    A(k)(s(j,jj,iii,jjjj)) = valu + A(k)(s(j,jj,iii,jjjj))
                  }
	        } // end jjjj 
	      } // end jj 
	    } // end j 

	    // for d^4u/dy^2dz^2 
            coefdyyzz = 6.0
            cyy = D.rderiv(2,ny,ii,hy)
            czz = D.rderiv(2,nz,iii,hz)
            for(jj <- 0 until ny )
	    {
              for(jjj <- 0 until nz )
              {
                valu = coefdyyzz*cyy(jj)*czz(jjj)
                if(jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1)
                {
	          A(k)(cs) = A(k)(cs) - valu*uana(xg(i),yg(jj),zg(jjj),tg(iiii))
                }
                else
                {
                  A(k)(s(i,jj,jjj,iiii)) = valu + A(k)(s(i,jj,jjj,iiii))
                }
	      } // end jjj 
	    } // end jj 

	    // for d^3u/dzdt^2 
            coefdztt  = 7.0
            cz = D.rderiv(1,nz,iii,hz)
            ctt = D.rderiv(2,nt,iiii,ht)
            for(jjj <- 0 until nz )
	    {
              for(jjjj <- 0 until nt )
              {
                valu = coefdztt*cz(jjj)*ctt(jjjj)
                if(jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1)
                {
	          A(k)(cs) = A(k)(cs) - valu*uana(xg(i),yg(ii),zg(jjj),tg(jjjj))
                }
                else
                { 
                  A(k)(s(i,ii,jjj,jjjj)) = valu + A(k)(s(i,ii,jjj,jjjj))
                }
	      } // end jjjj 
	    } // end jjj 

	    // for d^3u/dx^3 
            coefdxxx  = 8.0
            cxxx = D.rderiv(3,nx,i,hx)
            for(j <- 0 until nx )
            {
              valu = coefdxxx*cxxx(j)
              if(j==0 || j==nx-1)
              {
	        A(k)(cs) = A(k)(cs) - valu*uana(xg(j),yg(ii),zg(iii),tg(iiii))
              }
              else
              {
                A(k)(s(j,ii,iii,iiii)) = valu + A(k)(s(j,ii,iii,iiii))
              }
	    }

	    // for d^3u/dy^2dt 
            coefdyyt  = 9.0
            cyy = D.rderiv(2,ny,ii,hy)
            ct = D.rderiv(1,nt,iiii,ht)
            for(jj <- 0 until ny )
	    {
              for(jjjj <- 0 until nt )
              {
                valu = coefdyyt*cyy(jj)*ct(jjjj)
                if(jj==0 || jj==ny-1 || jjjj==0 || jjjj==nt-1)
                {
	          A(k)(cs) = A(k)(cs) - valu*uana(xg(i),yg(jj),zg(iii),tg(jjjj))
                }
                else
                {
                  A(k)(s(i,jj,iii,jjjj)) = valu +A(k)(s(i,jj,iii,jjjj))
                }
	      } // end jjjj 
	    } // end jj 

	    // for d^2u/dzdt 
            coefdzt   = 10.0
            cz = D.rderiv(1,nz,iii,hz)
            ct = D.rderiv(1,nt,iiii,ht)
            for(jjj <- 0 until nz )
	    {
              for(jjjj <- 0 until nt )
              {
                valu = coefdzt*cz(jjj)*ct(jjjj)
                if(jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1)
                {
	          A(k)(cs) = A(k)(cs) - valu*uana(xg(i),yg(ii),zg(jjj),tg(jjjj))
                }
                else
                {
                  A(k)(s(i,ii,jjj,jjjj)) = valu + A(k)(s(i,ii,jjj,jjjj))
                }
	      } // end jjjj 
	    } // end jjj 

	    // for du/dt 
            coefdt    = 11.0
            ct = D.rderiv(1,nt,iiii,ht)
            for(jjjj <- 0 until nt )
            {
              valu = coefdt*ct(jjjj)
              if(jjjj==0 || jjjj==nt-1)
              {
	        A(k)(cs) = A(k)(cs) - valu*uana(xg(i),yg(ii),zg(iii),tg(jjjj))
              }
              else
              {
                A(k)(s(i,ii,iii,jjjj)) = valu + A(k)(s(i,ii,iii,jjjj))
              }
	    }

	    // for u(x,y,z,t) 
            coefu     = 12.0
            A(k)(s(i,ii,iii,iiii)) = coefu + A(k)(s(i,ii,iii,iiii))

            // for f(x,y,z,t) RHS constant 
            A(k)(cs) = A(k)(cs) + f(xg(i),yg(ii),zg(iii),tg(iiii))

            // end terms 

  	  } // end iiii loop 
        } // end iii loop 
      } // end ii loop 
    } // end i loop 

    // printA() // lots of output 

    println("A computed stiffness matrix in "+(System.currentTimeMillis-Tmsec)+
            " milliseconds")
    println("")

    println("solve "+nxyzt+" equations in "+nxyzt+" unknowns ")

    U = S.solveb(A)

    println("equations solved at "+(System.currentTimeMillis-Tmsec)+
            " milliseconds")
    println("")
  
    println("check PDE against known solution")
    check_soln(Ua)
    println("check PDE against computed solution")
    check_soln(U)
    println("")

    println("          U computed, Ua analytic, error ")
    for(i <- 1 until nx-1 )
    {
      for(ii <- 1 until ny-1 )
      {
        for(iii <- 1 until nz-1 )
        {
          for(iiii <- 1 until nt-1 )
          {
            err = (U(s(i,ii,iii,iiii))-Ua(s(i,ii,iii,iiii))).abs
            maxerr = max(maxerr, err)
            avgerr = avgerr + err
            println("ug("+i+","+ii+","+iii+","+iiii+")="+
                    U(s(i,ii,iii,iiii))+",Ua="+
                    Ua(s(i,ii,iii,iiii))+", err="+err)
          } // iiii 
        } // iii 
      } // ii 
    } // i 
    //println("")
    println("   maxerr="+maxerr+", avgerr="+ avgerr/(nxyzt.toDouble) )
    Tmsec = System.currentTimeMillis-Tmsec
    println("total time "+Tmsec+" milliseconds")
    println("")
    println("end Pde44_eq.scala ")
  } // end main of Pde44_eq.scala


  def s(i:Int, ii:Int, iii:Int, iiii:Int):Int= // for x,y,z,t 
  {
    return (i-1)*(ny-2)*(nz-2)*(nt-2) + (ii-1)*(nz-2)*(nt-2) +
         (iii-1)*(nt-2) + (iiii-1)
  } // end s 

  def sk(k:Int, i:Int, ii:Int, iii:Int, iiii:Int):Int= // for x,y,z,t 
  {
    return k*(nxyzt+1) + (i-1)*(ny-2)*(nz-2)*(nt-2) + 
           (ii-1)*(nz-2)*(nt-2) + (iii-1)*(nt-2) + (iiii-1)
  } // end sk 

  def ss(i:Int, ii:Int, iii:Int, iiii:Int, j:Int, jj:Int, jjj:Int,
         jjjj:Int):Int =
  {
    return s(i,ii,iii,iiii)*(nxyzt+1) + s(j,jj,jjj,jjjj)
  } // end ss 

  def f(x:Double, y:Double, z:Double, t:Double):Double=
  {
    return 706.0 + 120.0*t*t + 140.0*y*y*z + 768.0*x + 180.0*z*z*t +
           200.0*y*y*z*t + 44.0*t*t*t + 110.0*y*y*z*z*t + 66.0*x*y +
           12.0*t*t*t*t + 24.0*z*z*z*z + 36.0*y*y*y*y + 48.0*x*x*x*x +
           60.0*y*y*z*z*t*t + 72.0*x*y*t + 84.0*x*z + 96.0*t 
  }

  def uana(x:Double, y:Double, z:Double, t:Double):Double=
                                      //analytic solution and boundaries
  {
    return t*t*t*t + 2.0*z*z*z*z + 3.0*y*y*y*y + 4.0*x*x*x*x +
           5.0*y*y*z*z*t*t + 6.0*x*y*t + 7.0*x*z + 8.0*t + 9.0
  } // end uana

  def printA()
  {
    var k:Int = 0
    var cs:Int = (nx-2)*(ny-2)*(nz-2)*(nt-2) // constant RHS column 
    for(i <- 1 until nx-1 )
    {
      for(ii <- 1 until ny-1 )
      {
        for(iii <- 1 until nz-1 )
        {
          for(iiii <- 1 until nt-1 )
          {
            k = s(i,ii,iii,iii) // row index 
            for(j <- 1 until nx-1 )
            {
              for(jj <- 1 until ny-1 )
              {
                for(jjj <- 1 until nz-1 )
                {
                  for(jjjj <- 1 until nt-1 )
                  {
                    println("A(row("+i+","+ii+","+iii+","+iiii+"),col("
                            +j+","+jj+","+jjj+","+jjjj+")) ="+
			    A(s(i,ii,iii,iiii))(s(j,jj,jjj,jjjj)))
	          } // jjjj 
	        } // jjj 
	      } // jj 
            } // j 
            println("RHS A("+k+","+cs+")="+A(k)(cs)) 
          } // iiii 
        } // iii 
      } // ii 
    } // i 
    println("")
  } // end printA 

  def evalu_derivative(xord:Int, yord:Int, zord:Int, tord:Int,
               i:Int, ii:Int, iii:Int, iiii:Int, u:Array[Double]):Double=
  {
    var cx=Array.ofDim[Double](nx)
    var cy=Array.ofDim[Double](ny)
    var cz=Array.ofDim[Double](nz)
    var ct=Array.ofDim[Double](nt)
    var p1=Array.ofDim[Double](nx)
    var p2=Array.ofDim[Double](ny)
    var p3=Array.ofDim[Double](nz)
    var discrete:Double = 0.0
    var x:Double = 0.0
    var y:Double = 0.0
    var z:Double = 0.0
    var t:Double = 0.0

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

    // four cases of one partial x, y, z, t to any order 
    if(xord!=0 && yord==0 && zord==0 && tord==0)
    {
      for(j <- 0 until nx )
      {  
        if(j==0 || j==nx-1)
        {
          discrete += cx(j)*uana(xg(j),y,z,t)
        }
        else
        {
          discrete += cx(j)*u(s(j,ii,iii,iiii))
        }
      }
    }
    else if(xord==0 && yord!=0 && zord==0 && tord==0)
      for(jj <- 0 until ny )
        if(jj==0 || jj==ny-1) discrete += cy(jj)*uana(x,yg(jj),z,t)
        else discrete += cy(jj)*u(s(i,jj,iii,iiii))
    else if(xord==0 && yord==0 && zord!=0 && tord==0)
      for(jjj <- 0 until nz )
        if(jjj==0 || jjj==nz-1) discrete += cz(jjj)*uana(x,y,zg(jjj),t)
        else discrete += cz(jjj)*u(s(i,ii,jjj,iiii))
    else if(xord==0 && yord==0 && zord==0 && tord!=0)
      for(jjjj <- 0 until nt )
        if(jjjj==0 || jjjj==nt-1) discrete += ct(jjjj)*uana(x,y,z,tg(jjjj))
        else discrete += ct(jjjj)*u(s(i,ii,iii,jjjj))

    // six cases of two partials xy, xz, xt, yz, yt, zt to any order 
    else if(xord!=0 && yord!=0 && zord==0 && tord==0)
    {
      for(j <- 0 until nx )
      {
        p1(j) = 0.0
        for(jj <- 0 until ny )
          if(j==0 || j==nx-1 || jj==0 || jj==ny-1)
               p1(j) += cy(jj)*uana(xg(j),yg(jj),zg(iii),tg(iiii))
          else   p1(j) += cy(jj)*u(s(j,jj,iii,iiii))
        discrete += cx(j)*p1(j)
      }
    }
    else if(xord!=0 && yord==0 && zord!=0 && tord==0)
    {
      for(j <- 0 until nx )
      {
        p1(j) = 0.0
        for(jjj <- 0 until nz )
          if(j==0 || j==nx-1 || jjj==0 || jjj==nz-1)
	       p1(j) += ct(jjj)*uana(xg(j),yg(ii),zg(jjj),tg(iiii))
          else p1(j) += cz(jjj)*u(s(j,ii,jjj,iiii))
        discrete += cx(j)*p1(j)
      }
    }
    else if(xord!=0 && yord==0 && zord==0 && tord!=0)
    {
      for(j <- 0 until nx )
      {
        p1(j) = 0.0
        for(jjjj <- 0 until nt )
          if(j==0 || j==nx-1 || jjjj==0 || jjjj==nt-1)
               p1(j) += ct(jjjj)*uana(xg(j),yg(ii),zg(iii),tg(jjjj))
          else p1(j) += ct(jjjj)*u(s(j,ii,iii,jjjj))
        discrete += cx(j)*p1(j)
      }
    }
    else if(xord==0 && yord!=0 && zord!=0 && tord==0)
    {
      for(jj <- 0 until ny )
      {
        p1(jj) = 0.0
        for(jjj <- 0 until nz )
          if(jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1)
               p1(jj) += cz(jjj)*uana(xg(i),yg(jj),zg(jjj),tg(iiii))
          else p1(jj) += cz(jjj)*u(s(i,jj,jjj,iiii))
        discrete += cy(jj)*p1(jj)
      }
    }
    else if(xord==0 && yord!=0 && zord==0 && tord!=0)
    {
      for(jj <- 0 until ny )
      {
        p1(jj) = 0.0
        for(jjjj <- 0 until nt )
          if(jj==0 || jj==ny-1 || jjjj==0 || jjjj==nt-1)
               p1(jj) += ct(jjjj)*uana(xg(i),yg(jj),zg(iii),tg(jjjj))
          else p1(jj) += ct(jjjj)*u(s(i,jj,iii,jjjj))
        discrete += cy(jj)*p1(jj)
      }
    }
    else if(xord==0 && yord==0 && zord!=0 && tord!=0)
    {
      for(jjj <- 0 until nz )
      {
        p1(jjj) = 0.0
        for(jjjj <- 0 until nt )
          if(jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1)
               p1(jjj) += ct(jjjj)*uana(xg(i),yg(ii),zg(jjj),tg(jjjj))
          else p1(jjj) += ct(jjjj)*u(s(i,ii,jjj,jjjj))
        discrete += cz(jjj)*p1(jjj)
      }
    }

    // four cases of three partials xyz, xyt, xzt, yzt to any order 
    else if(xord!=0 && yord!=0 && zord!=0 && tord==0)
    {
      for(j <- 0 until nx )
      {
        p1(j) = 0.0
        for(jj <- 0 until ny )
        {
          p2(jj) = 0.0
          for(jjj <- 0 until nz )
            if(j==0 || j==nx-1 ||
               jj==0 || jj==ny-1 ||
               jjj==0 || jjj==nz-1)
                 p2(jj) += cz(jjj)*uana(xg(j),yg(jj),zg(jjj),tg(iiii))
            else p2(jj) += cz(jjj)*u(s(j,jj,jjj,iiii))
          p1(j) += cy(jj)*p2(jj)
        }
        discrete += cx(j)*p1(j)
      }
    }
    else if(xord!=0 && yord!=0 && zord==0 && tord!=0)
    {
      for(j <- 0 until nx )
      {
        p1(j) = 0.0
        for(jj <- 0 until ny )
        {
          p2(jj) = 0.0
          for(jjjj <- 0 until nt )
            if(j==0 || j==nx-1 ||
               jj==0 || jj==ny-1 ||
               jjjj==0 || jjjj==nt-1)
                 p2(jj) += ct(jjjj)*uana(xg(j),yg(jj),zg(iii),tg(jjjj))
            else p2(jj) += ct(jjjj)*u(s(j,jj,iii,jjjj))
          p1(j) += cy(jj)*p2(jj)
        }
        discrete += cx(j)*p1(j)
      }
    }
    else if(xord!=0 && yord==0 && zord!=0 && tord!=0)
    {
      for(j <- 0 until nx )
      {
        p1(j) = 0.0
        for(jjj <- 0 until nz )
        {
          p2(jjj) = 0.0
          for(jjjj <- 0 until nt )
            if(j==0 || j==nx-1 ||
               jjj==0 || jjj==nz-1 ||
               jjjj==0 || jjjj==nt-1)
                 p2(jjj) += ct(jjjj)*uana(xg(j),yg(ii),zg(jjj),tg(jjjj))
            else p2(jjj) += ct(jjjj)*u(s(j,ii,jjj,jjjj))
          p1(j) += cz(jjj)*p2(jjj)
        }
        discrete += cx(j)*p1(j)
      }
    }
    else if(xord==0 && yord!=0 && zord!=0 && tord!=0)
    {
      for(jj <- 0 until ny )
      {
        p1(jj) = 0.0
        for(jjj <- 0 until nz )
        {
          p2(jjj) = 0.0
          for(jjjj <- 0 until nt )
            if(jj==0 || jj==ny-1 ||
               jjj==0 || jjj==nz-1 ||
               jjjj==0 || jjjj==nt-1)
                 p2(jjj) += ct(jjjj)*uana(xg(i),yg(jj),zg(jjj),tg(jjjj))
            else p2(jjj) += ct(jjjj)*u(s(i,jj,jjj,jjjj))
          p1(jj) += cz(jjj)*p2(jjj)
        }
        discrete += cy(jj)*p1(jj)
      }
    }
    // one case of four partials 
    else
    {
      for(j <- 0 until nx )
      {
        p1(j) = 0.0
        for(jj <- 0 until ny )
        {
          p2(jj) = 0.0
          for(jjj <- 0 until nz )
          {
            p3(jjj) = 0.0
            for(jjjj <- 0 until nt )
	      if(j==0 || j==nx-1 ||
                 jj==0 || jj==ny-1 ||
                 jjj==0 || jjj==nz-1 ||
                 jjjj==0 || jjjj==nt-1)
                 p3(jjj) += ct(jjjj)*uana(xg(j),yg(jj),zg(jjj),tg(jjjj))
	      else p3(jjj) += ct(jjjj)*u(s(j,jj,jjj,jjjj))
            p2(jj) += cz(jjj)*p3(jjj)
          } // end jjj 
          p1(j) += cy(jj)*p2(jj)
        } // end jj 
        discrete += cx(j)*p1(j)
      } // end j 
    } // end if 
    return discrete
  } // end evalu_deriv 

  def check_soln(u:Array[Double])
  {
    // check uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +
    //     4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +
    //     7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +
    //     10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) = F(x,y,z,t)
  
    var valu:Double = 0.0
    var err:Double = 0.0
    var smaxerr:Double = 0.0

    smaxerr = 0.0
    for(i <- 1 until nx-1 )
    {
      for(ii <- 1 until ny-1 )
      {
        for(iii <- 1 until nz-1 )
        {
          for(iiii <- 1 until nt-1 )
          {
            valu = 0.0
            // uxxxx(x,y,z,t) 
            valu = valu+evalu_derivative(4, 0, 0, 0, i, ii, iii, iiii, u)
	    // + 2 uyyyy(x,y,z,t) 
            valu = valu+2.0*evalu_derivative(0, 4, 0, 0, i, ii, iii, iiii, u)
	    // + 3 uzzzz(x,y,z,t) 
            valu = valu+3.0*evalu_derivative(0, 0, 4, 0, i, ii, iii, iiii, u)
            // + 4 utttt(x,y,z,t) 
            valu = valu+4.0*evalu_derivative(0, 0, 0, 4, i, ii, iii, iiii, u)
	    // + 5  uxyt(x,y,z,t) 
            valu = valu+5.0*evalu_derivative(1, 1, 0, 1, i, ii, iii, iiii, u)
	    // + 6 uyyzz(x,y,z,t) 
            valu = valu+6.0*evalu_derivative(0, 2, 2, 0, i, ii, iii, iiii, u)
            // + 7  uztt(x,y,z,t) 
            valu = valu+7.0*evalu_derivative(0, 0, 1, 2, i, ii, iii, iiii, u)
	    // + 8  uxxx(x,y,z,t) 
            valu = valu+8.0*evalu_derivative(3, 0, 0, 0, i, ii, iii, iiii, u)
	    // + 9  uyyt(x,y,z,t) 
            valu = valu+9.0*evalu_derivative(0, 2, 0, 1, i, ii, iii, iiii, u)
            // +10   uzt(x,y,z,t) 
            valu = valu+10.0*evalu_derivative(0, 0, 1, 1, i, ii, iii, iiii, u)
	    // +11    ut(x,y,z,t) 
            valu = valu+11.0*evalu_derivative(0, 0, 0, 1, i, ii, iii, iiii, u)
	    // +12     u(x,y,z,t) 
            valu = valu+12.0*uana(xg(i),yg(ii),zg(iii),tg(iiii))
            // -       F(x,y,z,t) should be zero 
	    err = (valu - f(xg(i),yg(ii),zg(iii),tg(iiii))).abs
            smaxerr = max(smaxerr, err)
          } // end iiii 
        } // end iii 
      } // end ii 
    } // end i 
    println("check_soln maxerr = ",smaxerr)
  } // end check_soln 
} // end object Pde44_eq