// pdeC44h_eq.java  homogeneous Biharmonic 4th order, 4 dimensions complex
// 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   u(x,y,z,t) = 0
//
// boundary conditions computed using Complex u(x,y,z,t)
// analytic solution is Complex u(x,y,z,t) = sin(x,y,z,t)
// uana(x,y,z,t) is used to computer boundary conditions
// and check accuracy of numeric solution
//
// 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 derivative order=4, npoints=5, point=i, term=0
// new Cnuderiv(order, npoints, point, xg[], cx[])
//
// d^4u/dx^4 using subscripts i,j,k,l 
// d^4u/dx^4 = sum(U(xg[i],yg[i],zg[k],tg[l])*cx[i]) i=0;i<nx
//
// subscripts in code use x[i], y[ii], z[iii], t[iiii]
// 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
// 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 matrix 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 side 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 Complex.java
// uses Csimeq.java
// uses Cinvert.java
// uses Cnuderiv.java

// example  Complex a, b, c;  c = a.add(b);  // does  c = a+b

class pdeC44h_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);

  Complex A[][] = new Complex[nxyzt][nxyzt+1]; // last column is RHS 
  Complex U[] = new Complex[nxyzt]; // computed solution
  Complex xg[] = new Complex[nx]; // grid values
  Complex yg[] = new Complex[ny];
  Complex zg[] = new Complex[nz];
  Complex tg[] = new Complex[nt];
  Complex Ua[] = new Complex[nxyzt]; // check solution
  Complex xmin = new Complex(0.0,0.1); // problem parameters 
  Complex xmax = new Complex(1.0,0.3);
  Complex ymin = new Complex(0.1,0.1);
  Complex ymax = new Complex(1.1,0.2);
  Complex zmin = new Complex(0.2,0.1);
  Complex zmax = new Complex(0.9,0.3);
  Complex tmin = new Complex(0.3,0.1);
  Complex tmax = new Complex(0.8,0.4);
  Complex hx, hy, hz, ht;
  Complex C2 = new Complex(2.0,0.0);

  pdeC44h_eq()
  {
    Complex x, y, z, t;
    int k, cs;
    // derivative point arrays, dynamic in this version 
    Complex cxx[] = new Complex[nx];
    Complex cyy[] = new Complex[ny];
    Complex czz[] = new Complex[nz];
    Complex ctt[] = new Complex[nt]; 
    Complex cxxxx[] = new Complex[nx];
    Complex cyyyy[] = new Complex[ny];
    Complex czzzz[] = new Complex[nz];
    Complex ctttt[] = new Complex[nt]; 
    Complex val;
    double  err, avgerr, maxerr;
    double time_start, time_now, time_last;

    System.out.println("pdeC44h_eq.java running ");
    System.out.println("Given uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) + ");
    System.out.println("    utttt(x,y,z,t) + 2 uxx(x,y,z,t) + 2 uyy(x,y,z,t) + ");
    System.out.println("    2 uzz(x,y,z,t) + 2 u(x,y,z,t) = 0");
    System.out.println("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries ");
    System.out.println("Analytic solution u(x,y,z,t) =  sin(x+y+z+t");
    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.subtract(xmin)).divide((double)(nx-1));
    for(int i=0; i<nx; i++)
    {
      xg[i] = xmin.add(hx.multiply((double)i));
      Ua[i] = uana(xg[i],ymin,zmin,tmin); // temp 
      System.out.println("i="+i+", Ua("+xg[i]+","+ymin+","+zmin+","+
                         tmin+")="+Ua[i]);
    } // i   
    System.out.println(" ");
    System.out.println("y grid and analytic solution at xmin,zmin,tmin ");
    hy = (ymax.subtract(ymin)).divide((double)(ny-1));
    for(int ii=0; ii<ny; ii++)
    {
      yg[ii] = ymin.add(hy.multiply((double)ii));
      Ua[ii] = uana(xmin, yg[ii],zmin,tmin); // temp 
      System.out.println("ii="+ii+", Ua("+xmin+","+yg[ii]+","+zmin+","+
                         tmin+")="+Ua[ii]);
    } // ii 
    System.out.println(" ");
    System.out.println("z grid and analytic solution ");
    hz = (zmax.subtract(zmin)).divide((double)(nz-1));
    for(int iii=0; iii<nz; iii++)
    {
      zg[iii] = zmin.add(hz.multiply((double)iii));
      Ua[iii] = uana(xmin, ymin, zg[iii],tmin); // temp 
    } // iii 
    System.out.println(" ");
    System.out.println("t grid and analytic solution ");
    ht = (tmax.subtract(tmin)).divide((double)(nt-1));
    for(int iiii=0; iiii<nt; iiii++)
    {
	tg[iiii] = tmin.add(ht.multiply((double)iiii));
      Ua[iiii] = uana(xmin, ymin, zmin, tg[iiii]); // temp 
    } // iiii 
    System.out.println(" ");

    System.out.println("check Cnuderiv(4,nx,0,xg,cxxxx");
    new Cnuderiv(4,nx,0,xg,cxxxx);
    for(int i=0; i<nx; i++) System.out.println("cxxxx["+i+"]="+cxxxx[i]);
    System.out.println("check Cnuderiv(4,nx,nx-1,xg,cxxxx");
    new Cnuderiv(4,nx,nx-1,xg,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 
    cs = nxyzt;
    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++)
          {
	    A[s(i,ii,iii,iiii)][cs] = new Complex(0.0,0.0);
            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)] = new Complex(0.0,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 = new Complex(0.0,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 
            new Cnuderiv(4,nx,i,xg,cxxxx);
            for(int j=0; j<nx; j++)
            {
              val = cxxxx[j];
              if(j==0 || j==nx-1)
		   A[k][cs] = A[k][cs].subtract(val.multiply(uana(xg[j],yg[ii],zg[iii],tg[iiii])));
              else A[k][s(j,ii,iii,iiii)] = A[k][s(j,ii,iii,iiii)].add(val);
	    }

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

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

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

	    // for 2 d^2u/dx^2 
            new Cnuderiv(2,nx,i,xg,cxx);
            for(int j=0; j<nx; j++)
            {
	      val = C2.multiply(cxx[j]);
              if(j==0 || j==nx-1)
		  A[k][cs] = A[k][cs].subtract(val.multiply(uana(xg[j],yg[ii],zg[iii],tg[iiii])));
              else A[k][s(j,ii,iii,iiii)] = A[k][s(j,ii,iii,iiii)].add(val);
	    }

	    // for 2 d^2u/dy^2 
            new Cnuderiv(2,ny,ii,yg,cyy);
            for(int jj=0; jj<ny; jj++)
            {
	      val = C2.multiply(cyy[jj]);
              if(jj==0 || jj==ny-1)
		  A[k][cs] =  A[k][cs].subtract(val.multiply(uana(xg[i],yg[jj],zg[iii],tg[iiii])));
              else A[k][s(i,jj,iii,iiii)] = A[k][s(i,jj,iii,iiii)].add(val);
	    }

	    // for 2 d^2u/dz^2 
            new Cnuderiv(2,nz,iii,zg,czz);
            for(int jjj=0; jjj<nz; jjj++)
            {
	      val = C2.multiply(czz[jjj]);
              if(jjj==0 || jjj==nz-1)
		  A[k][cs] = A[k][cs].subtract(val.multiply(uana(xg[i],yg[ii],zg[jjj],tg[iiii])));
              else A[k][s(i,ii,jjj,iiii)] = A[k][s(i,ii,jjj,iiii)].add(val);
	    }

	    // for u(x,y,z,t) 
            A[k][s(i,ii,iii,iiii)] = A[k][s(i,ii,iii,iiii)].add(C2);

            // 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 Csimeq(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 = (U[s(i,ii,iii,iiii)].subtract(Ua[s(i,ii,iii,iiii)])).abs();
            if(err>maxerr) maxerr = err;
            avgerr = avgerr + err;
            System.out.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 

    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="+maxerr+
		       ", avgerr="+(avgerr/(double)(nxyzt)));
    System.out.println(" ");
    System.out.println("end pdeC44h_eq.java");
  } // end pdeC44h_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, 0

  // problem boundaries (and analytic solution for checking)
  Complex uana(Complex x, Complex y, Complex z, Complex t)
  {
      return (x.add(y.add(z.add(t)))).sin();
  }

  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(Complex u[])
  {
    // 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 u(x,y,z,t) = 0
    Complex val;
    double 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 = new Complex(0.0,0.0);
            // uxxxx(x,y,z,t)
            val = val.add(eval_derivative(4, 0, 0, 0, i, ii, iii, iiii, u));
	    // + uyyyy(x,y,z,t)
            val = val.add(eval_derivative(0, 4, 0, 0, i, ii, iii, iiii, u));
	    // + uzzzz(x,y,z,t)
            val = val.add(eval_derivative(0, 0, 4, 0, i, ii, iii, iiii, u));
            // + utttt(x,y,z,t)
            val = val.add(eval_derivative(0, 0, 0, 4, i, ii, iii, iiii, u));
	    // + 2  uxx(x,y,z,t)
            val = val.add(C2.multiply(eval_derivative(2, 0, 0, 0, i, ii, iii, iiii, u)));
	    // + 2 uyy(x,y,z,t)
            val = val.add(C2.multiply(eval_derivative(0, 2, 0, 0, i, ii, iii, iiii, u)));
            // + 2 uzz(x,y,z,t)
            val = val.add(C2.multiply(eval_derivative(0, 0, 2, 0, i, ii, iii, iiii, u)));
	    // + 2 u(x,y,z,t)
            val = val.add(C2.multiply(uana(xg[i],yg[ii],zg[iii],tg[iiii])));
            //        F(x,y,z,t) should be zero
	    err = val.abs();
            if(err>smaxerr) smaxerr = err;
	  } // end iiii
	} // end iii
      } // end ii
    } // end i
    System.out.println("check_soln maxerr="+smaxerr);
  } // end check_soln

  Complex eval_derivative(int xord, int yord, int zord, int tord,
                         int i, int ii, int iii, int iiii, Complex u[])
  {
    Complex cx[] = new Complex[nx];
    Complex cy[] = new Complex[ny];
    Complex cz[] = new Complex[nz];
    Complex ct[] = new Complex[nt];
    int nn = Math.max(nx, Math.max(ny, Math.max(nz, nt)));
    Complex p1[] = new Complex[nn];
    Complex p2[] = new Complex[nn];
    Complex p3[] = new Complex[nn];
    Complex discrete;
    Complex x, y, z, t;

    new Cnuderiv(xord, nx, i, xg, cx);
    x = xg[i];
    new Cnuderiv(yord, ny, ii, yg, cy);
    y = yg[ii];
    new Cnuderiv(zord, nz, iii, zg, cz);
    z = zg[iii];
    new Cnuderiv(tord, nt, iiii, tg, ct);
    t = tg[iiii];
    discrete = new Complex(0.0,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 = discrete.add(cx[j].multiply(uana(xg[j],y,z,t)));
        else discrete = discrete.add(cx[j].multiply(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 = discrete.add(cy[jj].multiply(uana(x,yg[jj],z,t)));
        else discrete = discrete.add(cy[jj].multiply(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 = discrete.add(cz[jjj].multiply(uana(x,y,zg[jjj],t)));
        else discrete = discrete.add(cz[jjj].multiply(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 = discrete.add(ct[jjjj].multiply(uana(x,y,z,tg[jjjj])));
        else discrete = discrete.add(ct[jjjj].multiply(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] = new Complex(0.0,0.0);
        for(int jj=0; jj<ny; jj++)
	  if(j==0 || j==nx-1 || jj==0 || jj==ny-1)
	       p1[j] = p1[j].add(cy[jj].multiply(uana(xg[j],yg[jj],zg[iii],tg[iiii])));
          else p1[j] = p1[j].add(cy[jj].multiply(u[s(j,jj,iii,iiii)]));
        discrete = discrete.add(cx[j].multiply(p1[j]));
      } // end j
    }
    else if(xord!=0 && yord==0 && zord!=0 && tord==0)
    {
      for(int j=0; j<nx; j++)
      {
	p1[j] = new Complex(0.0,0.0);
        for(int jjj=0; jjj<nz; jjj++)
	  if(j==0 || j==nx-1 || jjj==0 || jjj==nz-1)
	       p1[j] = p1[j].add(ct[jjj].multiply(uana(xg[j],yg[ii],zg[jjj],tg[iiii])));
          else p1[j] = p1[j].add(cz[jjj].multiply(u[s(j,ii,jjj,iiii)]));
        discrete = discrete.add(cx[j].multiply(p1[j]));
      }
    }
    else if(xord!=0 && yord==0 && zord==0 && tord!=0)
    {
      for(int j=0; j<nx; j++)
      {
	p1[j] = new Complex(0.0,0.0);
        for(int jjjj=0; jjjj<nt; jjjj++)
	  if(j==0 || j==nx-1 || jjjj==0 || jjjj==nt-1)
	       p1[j] = p1[j].add(ct[jjjj].multiply(uana(xg[j],yg[ii],zg[iii],tg[jjjj])));
          else p1[j] = p1[j].add(ct[jjjj].multiply(u[s(j,ii,iii,jjjj)]));
        discrete = discrete.add(cx[j].multiply(p1[j]));
      }
    }
    else if(xord==0 && yord!=0 && zord!=0 && tord==0)
    {
      for(int jj=0; jj<ny; jj++)
      {
        p1[jj] = new Complex(0.0,0.0);
        for(int jjj=0; jjj<nz; jjj++)
	  if(jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1)
	       p1[jj] = p1[jj].add(cz[jjj].multiply(uana(xg[i],yg[jj],zg[jjj],tg[iiii])));
	  else p1[jj] = p1[jj].add(cz[jjj].multiply(u[s(i,jj,jjj,iiii)]));
        discrete = discrete.add(cy[jj].multiply(p1[jj]));
      }
    }
    else if(xord==0 && yord!=0 && zord==0 && tord!=0)
    {
      for(int jj=0; jj<ny; jj++)
      {
	p1[jj] = new Complex(0.0,0.0);
        for(int jjjj=0; jjjj<nt; jjjj++)
	  if(jj==0 || jj==ny-1 || jjjj==0 || jjjj==nt-1)
	       p1[jj] = p1[jj].add(ct[jjjj].multiply(uana(xg[i],yg[jj],zg[iii],tg[jjjj])));
          else p1[jj] = p1[jj].add(ct[jjjj].multiply(u[s(i,jj,iii,jjjj)]));
        discrete = discrete.add(cy[jj].multiply(p1[jj]));
      }
    }
    else if(xord==0 && yord==0 && zord!=0 && tord!=0)
    {
      for(int jjj=0; jjj<nz; jjj++)
      {
	p1[jjj] = new Complex(0.0,0.0);
        for(int jjjj=0; jjjj<nt; jjjj++)
	  if(jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1)
	       p1[jjj] = p1[jjj].add(ct[jjjj].multiply(uana(xg[i],yg[ii],zg[jjj],tg[jjjj])));
          else p1[jjj] = p1[jjj].add(ct[jjjj].multiply(u[s(i,ii,jjj,jjjj)]));
        discrete = discrete.add(cz[jjj].multiply(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] = new Complex(0.0,0.0);
        for(int jj=0; jj<ny; jj++)
	{
	  p2[jj] = new Complex(0.0,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] = p2[jj].add(cz[jjj].multiply(uana(xg[j],yg[jj],zg[jjj],tg[iiii])));
            else p2[jj] = p2[jj].add(cz[jjj].multiply(u[s(j,jj,jjj,iiii)]));
	  p1[j] = p1[j].add(cy[jj].multiply(p2[jj]));
	  }
        discrete = discrete.add(cx[j].multiply(p1[j]));
      }
    }
    else if(xord!=0 && yord!=0 && zord==0 && tord!=0)
    {
      for(int j=0; j<nx; j++)
      {
        p1[j] = new Complex(0.0,0.0);
        for(int jj=0; jj<ny; jj++)
	{
	  p2[jj] = new Complex(0.0,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] = p2[jj].add(ct[jjjj].multiply(uana(xg[j],yg[jj],zg[iii],tg[jjjj])));
            else p2[jj] = p2[jj].add(ct[jjjj].multiply(u[s(j,jj,iii,jjjj)]));
          p1[j] = p1[j].add(cy[jj].multiply(p2[jj]));
        }
        discrete = discrete.add(cx[j].multiply(p1[j]));
      }
    }
    else if(xord!=0 && yord==0 && zord!=0 && tord!=0)
    {
      for(int j=0; j<nx; j++)
      {
	p1[j] = new Complex(0.0,0.0);
        for(int jjj=0; jjj<nz; jjj++)
	{
	  p2[jjj] = new Complex(0.0,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] = p2[jjj].add(ct[jjjj].multiply(uana(xg[j],yg[ii],zg[jjj],tg[jjjj])));
            else p2[jjj] = p2[jjj].add(ct[jjjj].multiply(u[s(j,ii,jjj,jjjj)]));
          p1[j] = p1[j].add(cz[jjj].multiply(p2[jjj]));
	}
        discrete = discrete.add(cx[j].multiply(p1[j]));
      }
    }
    else if(xord==0 && yord!=0 && zord!=0 && tord!=0)
    {
      for(int jj=0; jj<ny; jj++)
      {
	p1[jj] = new Complex(0.0,0.0);
        for(int jjj=0; jjj<nz; jjj++)
	{
	  p2[jjj] = new Complex(0.0,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] = p2[jjj].add(ct[jjjj].multiply(uana(xg[i],yg[jj],zg[jjj],tg[jjjj])));
            else p2[jjj] = p2[jjj].add(ct[jjjj].multiply(u[s(i,jj,jjj,jjjj)]));
          p1[jj] = p1[jj].add(cz[jjj].multiply(p2[jjj]));
	} // end jjj
        discrete = discrete.add(cy[jj].multiply(p1[jj]));
      } // end jj
    }
    // one case of four partials
    else
    {
      for(int j=0; j<nx; j++)
      {
	p1[j] = new Complex(0.0,0.0);
        for(int jj=0; jj<ny; jj++)
	{
	  p2[jj] = new Complex(0.0,0.0);
          for(int jjj=0; jjj<nz; jjj++)
	  {
	    p3[jjj] = new Complex(0.0,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] = p3[jjj].add(ct[jjjj].multiply(uana(xg[j],yg[jj],zg[jjj],tg[jjjj])));
	      else p3[jjj] = p3[jjj].add(ct[jjjj].multiply(u[s(j,jj,jjj,jjjj)]));
            p2[jj] = p2[jj].add(cz[jjj].multiply(p3[jjj]));
	  } // end jjj
          p1[j] = p1[j].add(cy[jj].multiply(p2[jj]));
	} // end jj
        discrete = discrete.add(cx[j].multiply(p1[j]));
      } // end j
    } // end if
    return discrete;
  } // end eval_deriv

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