// pde44_eq.java  discretize, build linear equations, solve
// 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
 

// fourth order numerical difference order=4, npoints=5, term=0
// d^4u/dx^4 = (1/(b4[0]*hx^4))*(a4[0][0]*u(x,y,z,t) + a4[0][1]*u(x+hx,y,z,t) +
// a4[0][2]*u(x+2hx,y,z,t) + a4[0][3]*u(x+3hx,y,z,t) +a4[0][4]*u(x+4hx,y,z,t))
//
// d^4u/dx^4 using subscripts i,j,k,l for i a minimum subscript of x
// d^4u/dx^4 =(1/(b4[0]*hx^4))*(a4[0][0]*u(i,j,k,l) + a4[0][1]*u(i+1,j,k,l) +
// a4[0][2]*u(i+2,j,k,l) + a4[0][3]*u(i+3,j,k,l) +a4[0][4]*u(i+4,j,k,l))
//
// d^4u/dy^4 using subscripts i,j,k,l for j a minimum subscript of y, etc 
// d^4u/dy^4 =(1/(b4[0]*hy^4))*(a4[0][0]*u(i,j,k,l) + a4[0][1]*u(i,j+1,k,l) +
// a4[0][2]*u(i,j+2,k,l) + a4[0][3]*u(i,j+3,k,l) +a4[0][4]*u(i,j+4,k,l))
//
// d^4u/dx^4 using subscripts i,j,k,l for i-2 and i+2 allowed subscripts of x
// d^4u/dx^4 =(1/(b4[2]*hx^4))*(a4[2][0]*u(i-2,j,k,l) + a4[2][1]*u(i-1,j,k,l) +
// a4[2][2]*u(i,j,k,l) + a4[2][3]*u(i+1,j,k,l) +a4[2][4]*u(i+2,j,k,l))
//
// d^4u/dx^4 using subscripts i,j,k,l for i a maximum subscript of x
// d^4u/dx^4 =(1/(b4[4]*hx^4))*(a4[4][0]*u(i-4,j,k,l) + a4[4][1]*u(i-3,j,k,l) +
// a4[4][2]*u(i-2,j,k,l) + a4[4][3]*u(i-1,j,k,l) +a4[4][4]*u(i,j,k,l))
//
// 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
 
// uses simeq.java
// uses deriv.java
import java.text.*;

class pde44_eq
{
  final int nx = 6; // problem parameters
  final int ny = 6;
  final int nz = 6;
  final int nt = 6;
  final int nxyzt = (nx-2)*(ny-2)*(nz-2)*(nt-2);

  double A[][] = new double[nxyzt][nxyzt+1]; // last column is RHS 
  double U[] = new double[nxyzt];
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double tg[] = new double[nt];
  double Ua[] = new double[nxyzt];
  double xmin = 0.0; // problem parameters 
  double xmax = 1.0;
  double ymin = 0.0;
  double ymax = 1.0;
  double zmin = 0.0;
  double zmax = 1.0;
  double tmin = 0.0;
  double tmax = 1.0;
  double hx, hy, hz, ht;

  DecimalFormat f6 = new DecimalFormat("0.000");
  DecimalFormat fe = new DecimalFormat("0.000E00");

  pde44_eq()
  {
    double x, y, z, t;
    int k, cs;
    // coeficients of terms in first equation 
    double coefdxxxx; // computed below if not constant 
    double coefdyyyy;
    double coefdzzzz;
    double coefdtttt;
    double coefdxyt;
    double coefdyyzz;
    double coefdztt;
    double coefdxxx;
    double coefdyyt;
    double coefdzt;
    double coefdt;
    double coefu;
    // derivative point arrays, dynamic in this version 
    double cx[] = new double[nx];
    double cy[] = new double[ny];
    double cz[] = new double[nz];
    double ct[] = new double[nt]; 
    double cxx[] = new double[nx];
    double cyy[] = new double[ny];
    double czz[] = new double[nz];
    double ctt[] = new double[nt]; 
    double cxxx[] = new double[nx];
    double cyyy[] = new double[ny];
    double czzz[] = new double[nz];
    double cttt[] = new double[nt]; 
    double cxxxx[] = new double[nx];
    double cyyyy[] = new double[ny];
    double czzzz[] = new double[nz];
    double ctttt[] = new double[nt]; 
    double val, err, avgerr, maxerr;
    double time_start;
    double time_now;
    double time_last;

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

    System.out.println("xmin="+xmin+", xmax="+xmax+", ymin="+ymin+", ymax="+ymax);
    System.out.println("zmin="+zmin+", zmax="+zmax+", tmin="+tmin+", tmax="+tmax);
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz+", nt="+nt);
    time_start = (double)System.currentTimeMillis()/1000.0;
    time_last = time_start;

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

    System.out.println("check rderiv(4,nx,0,hx,cxxxx");
    new rderiv(4,nx,0,hx,cxxxx);
    for(int i=0; i<nx; i++) System.out.println("cxxxx["+i+"]="+cxxxx[i]);
    System.out.println("check rderiv(4,nx,nx-1,hx,cxxxx");
    new rderiv(4,nx,nx-1,hx,cxxxx);
    for(int i=0; i<nx; i++) System.out.println("cxxxx["+i+"]="+cxxxx[i]);
    System.out.println(" ");

    System.out.println("put solution in Ua vector ");
    // put solution in Ua vector 
    for(int i=1; i<nx-1; i++)
    {
      x = xg[i];
      for(int ii=1; ii<ny-1; ii++)
      {
        y = yg[ii];
        for(int iii=1; iii<nz-1; iii++)
        {
          z = zg[iii];
          for(int iiii=1; iiii<nt-1; iiii++)
          {
            t = tg[iiii];
            Ua[s(i,ii,iii,iiii)] = uana(x,y,z,t);
          } // iiii 
        } // iii 
      } // ii 
    } // i 
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("time for previous section = "+(time_now-time_last)+
                       " seconds ");
    time_last = time_now;
    System.out.println(" ");

    System.out.println("initialize A ");
    // initialize A 
    for(int i=1; i<nx-1; i++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        for(int iii=1; iii<nz-1; iii++)
        {
          for(int iiii=1; iiii<nt-1; iiii++)
          {
            for(int j=1; j<nx-1; j++)
            {
              for(int jj=1; jj<ny-1; jj++)
              {
                for(int jjj=1; jjj<nz-1; jjj++)
                {
                  for(int jjjj=1; jjjj<nt-1; jjjj++)
                  {
		    A[s(i,ii,iii,iiii)][s(j,jj,jjj,jjjj)] = 0.0;
	          } // jjjj 
	        } // jjj 
	      } // jj 
            } // j 
          } // iiii 
        } // iii 
      } // ii 
    } // i 

    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("time for previous section = "+(time_now-time_last)+
                       " seconds ");
    time_last = time_now;
    System.out.println(" ");

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

            // for each term in first equation 

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

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

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

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

	    // for d^3u/dxdydt 
            coefdxyt  = 5.0;
            new rderiv(1,nx,i,hx,cx);
            new rderiv(1,ny,ii,hy,cy);
            new rderiv(1,nt,iiii,ht,ct);
            for(int j=0; j<nx; j++)
	    {
              for(int jj=0; jj<ny; jj++)
	      {
                for(int jjjj=0; jjjj<nt; jjjj++)
                {
                  val = 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] -= val*uana(xg[j],yg[jj],zg[iii],tg[jjjj]);
                  else A[k][s(j,jj,iii,jjjj)] += val;
	        } // end jjjj 
	      } // end jj 
	    } // end j 

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

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

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

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

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

	    // for du/dt 
            coefdt    = 11.0;
            new rderiv(1,nt,iiii,ht,ct);
            for(int jjjj=0; jjjj<nt; jjjj++)
            {
              val = coefdt*ct[jjjj];
              if(jjjj==0 || jjjj==nt-1)
	           A[k][cs] -= val*uana(xg[i],yg[ii],zg[iii],tg[jjjj]);
              else A[k][s(i,ii,iii,jjjj)] += val;
	    }

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

            // for f(x,y,z,t) RHS constant 
            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();

    System.out.println("A computed stiffness matrix ");
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("time for previous section = "+(time_now-time_last)+
                       " seconds ");
    time_last = time_now;
    System.out.println(" ");


    System.out.println("solve "+nxyzt+" equations in "+nxyzt+" unknowns ");

    new simeq(nxyzt, A, U);

    System.out.println("equations solved ");
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("time for previous section = "+(time_now-time_last)+
                       " seconds ");
    time_last = time_now;
    System.out.println(" ");

    System.out.println("check PDE against known solution");
    check_soln(Ua);
    System.out.println("check PDE against computed solution");
    check_soln(U);
    System.out.println(" ");

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("          U computed, Ua analytic, error ");
    for(int i=1; i<nx-1; i++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        for(int iii=1; iii<nz-1; iii++)
        {
          for(int iiii=1; iiii<nt-1; iiii++)
          {
            err = Math.abs(U[s(i,ii,iii,iiii)]-Ua[s(i,ii,iii,iiii)]);
            if(err>maxerr) maxerr = err;
            avgerr = avgerr + err;
            System.out.println("ug["+i+","+ii+","+iii+","+iiii+"]="+
			       f6.format(U[s(i,ii,iii,iiii)])+", Ua="+
                               f6.format(Ua[s(i,ii,iii,iiii)])+", err="+
                               fe.format(err));
          } // iiii 
        } // iii 
      } // ii 
    } // i 

    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("total CPU time  = "+(time_now-time_start)+" seconds"); 
    System.out.println(" ");


    System.out.println("                            maxerr="+fe.format(maxerr)+
		       ", avgerr="+fe.format((avgerr/(double)(nx*ny*nz*nt))));
    System.out.println(" ");
    System.out.println("end pde44_eq.java");
  } // end pde44_eq

  int s(int i, int ii, int iii, int iiii) // 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 

  int sk(int k, int i, int ii, int iii, int iiii) // 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 

  int ss(int i, int ii, int iii, int iiii, int j, int jj, int jjj, int jjjj)
  {
    return s(i,ii,iii,iiii)*(nxyzt+1) + s(j,jj,jjj,jjjj);
  } // end ss 

  // problem RHS
  double f(double x, double y, double z, double t)
  {
    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; 
  }

  // problem boundaries (and analytic solution for checking)
  double uana(double x, double y, double z, double t)
  {
    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;
  }

  void printA()
  {
    int cs, k;
    cs = (nx-2)*(ny-2)*(nz-2)*(nt-2); // constant RHS column 
    for(int i=1; i<nx-1; i++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        for(int iii=1; iii<nz-1; iii++)
        {
          for(int iiii=1; iiii<nt-1; iiii++)
          {
            k = s(i,ii,iii,iiii); // row index 
            for(int j=1; j<nx-1; j++)
            {
              for(int jj=1; jj<ny-1; jj++)
              {
                for(int jjj=1; jjj<nz-1; jjj++)
                {
                  for(int jjjj=1; jjjj<nt-1; jjjj++)
                  {
                    System.out.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 
            System.out.println("RHS A["+k+","+cs+"]="+A[k][cs]);
            // equivalent A[sk(k,i,ii,iii,iiii+1)] 
          } // iiii 
        } // iii 
      } // ii 
    } // i 
    System.out.println(" ");
  } // end printA 

  void check_soln(double u[])
  {
    // 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)
    double val, err, smaxerr;

    smaxerr = 0.0;
    for(int i=1; i<nx-1; i++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        for(int iii=1; iii<nz-1; iii++)
        {
          for(int iiii=1; iiii<nt-1; iiii++)
          {
            val = 0.0;
            // uxxxx(x,y,z,t)
            val += eval_derivative(4, 0, 0, 0, i, ii, iii, iiii, u);
	    // + 2 uyyyy(x,y,z,t)
            val += 2.0*eval_derivative(0, 4, 0, 0, i, ii, iii, iiii, u);
	    // + 3 uzzzz(x,y,z,t)
            val += 3.0*eval_derivative(0, 0, 4, 0, i, ii, iii, iiii, u);
            // + 4 utttt(x,y,z,t)
            val += 4.0*eval_derivative(0, 0, 0, 4, i, ii, iii, iiii, u);
	    // + 5  uxyt(x,y,z,t)
            val += 5.0*eval_derivative(1, 1, 0, 1, i, ii, iii, iiii, u);
	    // + 6 uyyzz(x,y,z,t)
            val += 6.0*eval_derivative(0, 2, 2, 0, i, ii, iii, iiii, u);
            // + 7  uztt(x,y,z,t)
            val += 7.0*eval_derivative(0, 0, 1, 2, i, ii, iii, iiii, u);
	    // + 8  uxxx(x,y,z,t)
            val += 8.0*eval_derivative(3, 0, 0, 0, i, ii, iii, iiii, u);
	    // + 9  uyyt(x,y,z,t)
            val += 9.0*eval_derivative(0, 2, 0, 1, i, ii, iii, iiii, u);
            // +10   uzt(x,y,z,t)
            val += 10.0*eval_derivative(0, 0, 1, 1, i, ii, iii, iiii, u);
	    // +11    ut(x,y,z,t)
            val += 11.0*eval_derivative(0, 0, 0, 1, i, ii, iii, iiii, u);
	    // +12     u(x,y,z,t)
            val += 12.0*uana(xg[i],yg[ii],zg[iii],tg[iiii]);
            // -       F(x,y,z,t) should be zero
	    err = Math.abs(val - f(xg[i],yg[ii],zg[iii],tg[iiii]));
            if(err>smaxerr) smaxerr = err;
	  } // end iiii
	} // end iii
      } // end ii
    } // end i
    System.out.println("check_soln maxerr="+smaxerr);
  } // end check_soln
 
  double eval_derivative(int xord, int yord, int zord, int tord,
                         int i, int ii, int iii, int iiii, double u[])
  {
    double cx[] = new double[nx];
    double cy[] = new double[ny];
    double cz[] = new double[nz];
    double ct[] = new double[nt];
    int nn = Math.max(nx, Math.max(ny, Math.max(nz, nt)));
    double p1[] = new double[nn];
    double p2[] = new double[nn];
    double p3[] = new double[nn];
    double discrete;
    double x, y, z, t;

    new rderiv(xord, nx, i, hx, cx);
    x = xg[i];
    new rderiv(yord, ny, ii, hy, cy);
    y = yg[ii];
    new rderiv(zord, nz, iii, hz, cz);
    z = zg[iii];
    new rderiv(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!=0 && yord==0 && zord==0 && tord==0)
      for(int j=0; j<nx; j++)
	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(int jj=0; jj<ny; jj++)
	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(int jjj=0; jjj<nz; jjj++)
	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(int jjjj=0; jjjj<nt; jjjj++)
	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(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jj=0; jj<ny; jj++)
	  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];
      } // end j
    }
    else if(xord!=0 && yord==0 && zord!=0 && tord==0)
    {
      for(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jjj=0; jjj<nz; jjj++)
	  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(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jjjj=0; jjjj<nt; jjjj++)
	  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(int jj=0; jj<ny; jj++)
      {
        p1[jj] = 0.0;
        for(int jjj=0; jjj<nz; jjj++)
	  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(int jj=0; jj<ny; jj++)
      {
        p1[jj] = 0.0;
        for(int jjjj=0; jjjj<nt; jjjj++)
	  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(int jjj=0; jjj<nz; jjj++)
      {
        p1[jjj] = 0.0;
        for(int jjjj=0; jjjj<nt; jjjj++)
	  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(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jj=0; jj<ny; jj++)
	{
          p2[jj] = 0.0;
          for(int jjj=0; jjj<nz; jjj++)
	    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(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jj=0; jj<ny; jj++)
	{
          p2[jj] = 0.0;
          for(int jjjj=0; jjjj<nt; jjjj++)
	    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(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jjj=0; jjj<nz; jjj++)
	{
          p2[jjj] = 0.0;
          for(int jjjj=0; jjjj<nt; jjjj++)
	    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(int jj=0; jj<ny; jj++)
      {
        p1[jj] = 0.0;
        for(int jjj=0; jjj<nz; jjj++)
	{
          p2[jjj] = 0.0;
          for(int jjjj=0; jjjj<nt; jjjj++)
	    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];
	} // end jjj
        discrete += cy[jj]*p1[jj];
      } // end jj
    }
    // one case of four partials
    else
    {
      for(int j=0; j<nx; j++)
      {
        p1[j] = 0.0;
        for(int jj=0; jj<ny; jj++)
	{
          p2[jj] = 0.0;
          for(int jjj=0; jjj<nz; jjj++)
	  {
            p3[jjj] = 0.0;
            for(int jjjj=0; jjjj<nt; jjjj++)
		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 eval_deriv

  public static void main (String[] args)
  {
    new pde44_eq();
  }
} // end class pde44_eq