// pde24js_eq.java  2D uniform grid boundary value PDE  based on pde24js_eq.js
//                known solution to test method
// 
// solve equation in uana24.js
// uana24(x,y,z,t), f(x,y,z,t)
//
// set up and solve system of (nx-2)*(ny-2)*(nz-2)*(nt-2) linear equations
//
// uses new nuderiv.java
// uses simeq.java
// uses uana24.js  uana()  f()  coefs()  define problem 

class pde24js_eq
{
pde24js_eq()
{
  double start, now;
  boolean debug = false;
  double xg[] = new double[10];   // data points
  double yg[] = new double[10];
  double zg[] = new double[10];
  double tg[] = new double[10];
  int neqn; 
  int cs;
  double us[] = new double[10000];   // solution being computed 
  double ut[][] = new double[10000][10001] ;   // [neqn][neqn+1];  equations incl RHS 
  double coef[] = new double[29];
  double ccxx, ccyy, cczz, cctt, ccxy, ccxz, ccxt, ccyz, ccyt, cczt;
  double ccx, ccy, ccz, cct, ccc;
  int nx, ny, nz, nt;
  double xmin, xmax, hx;
  double ymin, ymax, hy;
  double zmin, zmax, hz;
  double tmin, tmax, ht;
  double error;         // sum of absolute exact-computed         
  double avgerror;      // average error                          
  double maxerror;      // maximum cell error                     

  System.out.print("pde24js_eq.java running \n");
  System.out.print("differential equation to solve, from uana24.java \n");
  System.out.print("given f(x,y,z), uana(x,y,x) \n");
  System.out.print("uniform grid on rectangle xmin,ymin,zmin to xmax,ymax,zmax\n");
  System.out.print(" \n");

  ccxx = 1.0;
  ccyy = 2.0;
  cczz = 3.0;
  cctt = 4.0;
  ccxy = 5.0;
  ccxz = 6.0;
  ccxt = 7.0;
  ccyz = 8.0;
  ccyt = 9.0;
  cczt = 10.0;
  ccx  = 11.0;
  ccy  = 12.0;
  ccz  = 13.0;
  cct  = 14.0;
  ccc  = 15.0;
  xmin = 0.0;
  xmax = 1.0;
  nx   = 5;
  ymin = 0.1;
  ymax = 1.1;
  ny   = 6;
  zmin = 0.2;
  zmax = 1.2;
  nz   = 7;
  tmin = 0.3;
  tmax = 1.3;
  nt   = 8;
  System.out.print("ccxx="+ccxx+", ccyy="+ccyy+", cczz="+cczz+" \n");
  System.out.print("cctt="+cctt+", ccxy="+ccxy+", ccxz="+ccxz+" \n");
  System.out.print("ccxt="+ccxt+", ccyz="+ccyz+", ccyt="+ccyt+" \n");
  System.out.print("cczt="+cczt+" \n");
  System.out.print("ccx="+ccx+", ccy="+ccy+", ccz="+ccz+" \n");
  System.out.print("cct="+cct+", ccc="+ccc+" \n");
  System.out.print("xmin="+xmin+", xmax="+xmax+", nx="+nx+" \n");
  System.out.print("ymin="+ymin+", ymax="+ymax+", ny="+ny+" \n");
  System.out.print("zmin="+zmin+", zmax="+zmax+", nz="+nz+" \n");
  System.out.print("tmin="+tmin+", tmax="+tmax+", nt="+nt+" \n");
  start = System.currentTimeMillis()/1000.0;

  neqn = (nx-2)*(ny-2)*(nz-2)*(nt-2); // internal, non boundary cells 
  System.out.print(neqn+" equations will be set up and solved \n");
  cs = neqn;
  hx = (xmax-xmin)/(nx-1);
  System.out.print("hx="+hx+" \n"); 

  for(int i=0; i<nx; i++)
  {
    xg[i] = xmin+i*hx;
      System.out.print("xg["+i+"]="+xg[i]+" \n");
  }

  hy = (ymax-ymin)/(ny-1);
  System.out.print("hy="+hy+" \n"); 
  for(int j=0; j<ny; j++)
  {
    yg[j] = ymin+j*hy;
    System.out.print("yg("+j+")="+yg[j]+" \n");
  }

  hz = (zmax-zmin)/(nz-1);
  System.out.print("hz="+hz+" \n"); 
  for(int k=0; k<nz; k++)
  {
    zg[k] = zmin+k*hz;
    System.out.print("zg("+k+")="+zg[k]+" \n");
  }

  ht = (tmax-tmin)/(nt-1);
  System.out.print("ht="+ht+" \n"); 
  for(int i4=0; i4<nt; i4++)
  {
    tg[i4] = tmin+i4*ht;
    System.out.print("tg("+i4+")="+tg[i4]+" \n");
  }

  System.out.print("test uana(xmax,ymax,zmax,tmax)="+
                   uana(xmax,ymax,zmax,tmax)+" \n");
  System.out.print("test f(xmax,ymax,zmax,tmax)="+
                   f(xmax,ymax,zmax,tmax)+" \n");

  for(int i=0; i<neqn; i++)
  {
    for(int j=0; j<=neqn; j++)
	ut[i][j] = 0.0;
  }
  System.out.print("ut initialized \n");
  init_matrix(nx, ny, nz, nt, ut, xg, yg, zg, tg,
              ccxx, ccyy, cczz, cctt, ccxy, ccxz, ccxt,
              ccyz, ccyt, cczt, ccx, ccy, ccz, cct, ccc);
  System.out.print("matrix initialized \n");
  if(debug) printB(nx, ny, nz, nt, ut); // for debugging
  check_rows(nx, ny, nz, nt, ut, xg, yg, zg, tg); 

  System.out.print("running simeq \n");
  new simeq(neqn, ut, us);
  System.out.print("back from simeq \n");

  check_soln(nx, ny, nz, nt, us, xg, yg, zg, tg,
             ccxx, ccyy, cczz, cctt, ccxy, ccxz, ccxt,
             ccyz, ccyt, cczt, ccx, ccy, ccz, cct, ccc);
  printS(nx, ny, nz, nt, xg, yg, zg, tg, us);

  check(nx, ny, nz, nt, xg, yg, zg, tg, us);

  now = System.currentTimeMillis()/1000.0;
  System.out.println("total time "+(now-start)+" seconds");
  System.out.print("finished pde24js_eq.java \n");
} // end pde24js_eq

int S(int i, int j, int k, int i4, int ny, int nz, int nt)
                            // zero based index includes boundary 
{                           // this is row of ut[] 
  return (i-1)*(ny-2)*(nz-2)*(nt-2)+(j-1)*(nz-2)*(nt-2)+(nt-2)*(k-1)+(i4-1); 
         // index into us and ut, RHS, no boundary 
} // end S 

void check(int nx, int ny, int nz, int nt, double xg[], double yg[],
           double zg[], double tg[], double us[]) 
              // typically not known, instrumentation 
{
  double uexact, err, error, maxerror;
  
  error = 0.0;
  maxerror = 0.0;
  for(int i=1; i<nx-1; i++)
  {
    for(int j=1; j<ny-1; j++)
    { 
      for(int k=1; k<nz-1; k++)
      {
        for(int i4=1; i4<nt-1; i4++)
        {
          uexact = uana(xg[i], yg[j], zg[k], tg[i4]);
          err = Math.abs(us[S(i,j,k,i4,ny,nz,nt)]-uexact);
          error = error + err;
          if(err>maxerror) maxerror=err;
        } // i4
      } // k
    } // j
  } // i
  System.out.print("check error="+error+", maxerror="+maxerror+" \n");
} // end check 

void check_rows(int nx, int ny, int nz, int nt,
                double ut[][], double xg[], double yg[],
                double zg[], double tg[]) // only for(checking 
{	// row is in ut solution coordinates 
  int row, col;
  double sum, sum1;
  int cs = (nx-2)*(ny-2)*(nz-2)*(nt-2);

  System.out.print("in check_rows \n");
  for(int i=1; i<nx-1; i++)
  {
    for(int j=1; j<ny-1; j++)
    {
      for(int k=1; k<nz-1; k++)
      {
        for(int i4=1; i4<nt-1; i4++)
        {
	  row = S(i,j,k,i4,ny,nz,nt);
          sum1 = 0.0;
          for(int ii=1; ii<nx-1; ii++)
          {
            for(int jj=1; jj<ny-1; jj++)
	    {
              for(int kk=1; kk<nz-1; kk++)
              {
                for(int j4=1; j4<nt-1; j4++)
                {
                  col = S(ii,jj,kk,j4,ny,nz,nt);
                  sum1 = sum1 + ut[row][col]*uana(xg[ii],yg[jj],zg[kk],tg[j4]);
                            // u is check solution 
                } // j4 
              } // kk
            } // jj
          } // ii
          sum = sum1 - ut[row][cs];
          if(Math.abs(sum)>0.001)
          {
            System.out.print("row i="+i+", j="+j+", k="+k+", i4="+i4+" is bad err="+sum+" \n");
            System.out.print("ut[row][cs]="+ut[row][cs]+", sum1="+sum1+" \n"); 
          } // end if
        } // i4
      } // k
    } // j 
  } //i 
  System.out.print("finished check_rows \n");
} // end check_rows 

void printS(int nx, int ny, int nz, int nt, double xg[], double yg[],
            double zg[], double tg[], double us[]) // typically not known 
{
  double error = 0.0;
  double max_error = 0.0;
  double avg_error = 0.0;
  int neqn = (nx-2)*(ny-2)*(nz-2)*(nt-2);

  System.out.print("exact solution u, computed us, error \n");
  for(int i=1; i<nx-1; i++)
  {
    for(int j=1; j<ny-1; j++)
    {
      for(int k=1; k<nz-1; k++)
      {
        for(int i4=1; i4<nt-1; i4++)
        {
          error = uana(xg[i],yg[j],zg[k],tg[i4]) - us[S(i,j,k,i4,ny,nz,nt)];
          avg_error = avg_error + Math.abs(error);
          if(Math.abs(error)>max_error) max_error = Math.abs(error);

          System.out.print("xg["+i+"]="+xg[i]+", yg["+j+"]="+yg[j]+
                         ", zg["+k+"]="+zg[k]+", tg["+i4+"]="+tg[i4]+" \n");
	  System.out.print("u="+uana(xg[i],yg[j],zg[k],tg[i4])+
                         ", us="+us[S(i,j,k,i4,ny,nz,nt)]+
                         ", err="+error+" \n");
        } // i4
      } // k
    } // j
  } // i
  System.out.print("max_error="+max_error+", avg_error="+avg_error/neqn+" \n");
  System.out.print(" \n");
} // end print

void printB(int nx, int ny, int nz, int nt, double ut[][])
{
  System.out.print("matrix B includes RHS \n");
  int cs = (nx-2)*(ny-2)*(nz-2)*(nt-2);
  for(int i=1; i<nx-1; i++)
  {
    for(int j=1; j<ny-1; j++)
    {
      for(int k=1; k<nz-1; k++)
      {
        for(int i4=1; i4<nt-1; i4++)
        {
          for(int ii=1; ii<nx-1; ii++)
          {
            for(int jj=1; jj<ny-1; jj++)
            {
              for(int kk=1; kk<nz-1; kk++)
              {
                for(int j4=1; j4<nt-1; j4++)
                {
                  System.out.print("i="+i+", j="+j+", k="+k+", i4="+i4+
                                 ", ii="+ii+", jj="+jj+", kk="+kk+
                                 ", j4="+j4+", ut="+
	                    ut[S(i,j,k,i4,ny,nz,nt)][S(ii,jj,kk,j4,ny,nz,nt)]+
                            " \n");
                } // j4
              } // kk
	    } // jj 
          } // ii 
          System.out.print("RHS=           "+ut[S(i,j,k,i4,ny,nz,nt)][cs]+" \n");
        } // i4
      } // k
    } // j 
  } // i 
  System.out.print(" \n");
  System.out.print(" \n");
} //  end printB 

void init_matrix(int nx, int ny, int nz, int nt, double ut[][],
                 double xg[], double yg[], double zg[], double tg[],
                 double ccxx, double ccyy, double cczz, double cctt,
                 double ccxy, double ccxz, double ccxt,
                 double ccyz, double ccyt, double cczt,
                 double ccx,  double ccy,  double ccz, double cct, double ccc)
{
  double cxx[] = new double[10];  // new nuderiv coefficients 
  double cyy[] = new double[10];
  double czz[] = new double[10];
  double ctt[] = new double[10];
  double cxy[] = new double[10];
  double cxz[] = new double[10];
  double cxt[] = new double[10];
  double cyz[] = new double[10];
  double cyt[] = new double[10];
  double czt[] = new double[10];
  double cx[] = new double[10]; 
  double cy[] = new double[10];
  double cz[] = new double[10];
  double ct[] = new double[10];
  int row;
  int cs = (nx-2)*(ny-2)*(nz-2)*(nt-2);

  System.out.print("init_matrix running nx="+nx+", ny="+ny+", nz="+nz+" \n");
  System.out.print("nt="+nt+", cs="+cs+" \n");
  System.out.print("S(1,1,1,1,ny,nz,nt)="+S(1,1,1,1,ny,nz,nt)+" \n");
  System.out.print("S(nx-2,ny-2,nz-2,nt-2,ny,nz,nt)="+
                  S(nx-2,ny-2,nz-2,nt-2,ny,nz,nt)+" \n");

  // cs is RHS, constant column subscript 
  // zero internal cells 
  for(int i=1; i<nx-1; i++)
  {
    for(int j=1; j<ny-1; j++)
    {
      for(int k=1; k<nz-1; k++)
      {
        for(int i4=1; i4<nt-1; i4++)
        {
          row = S(i,j,k,i4,ny,nz,nt);
          // System.out.print("row S("+i+","+j+","+k+","+i4+","+
          //                ny+","+nz+","+nt+")="+
          //                S(i,j,k,ny,nz)+" \n");
          for(int ii=1; ii<nx-1; ii++)
          {
            for(int jj=1; jj<ny-1; jj++)
            {
              for(int kk=1; kk<nz-1; kk++)
              {
                for(int j4=1; j4<nt-1; j4++)
                {
	          ut[row][S(ii,jj,kk,j4,ny,nz,nt)] = 0.0;
                } // j4
              } // kk
            } //jj
          } // ii
          ut[row][cs] = 0.0;
        } // i4
      } // k
    } // j
  } // i

  System.out.print("internal cells zeroed \n");
  for(int i=1; i<nx-1; i++) // inside boundary 
  {
    for(int j=1; j<ny-1; j++)
    {
      for(int k=1; k<nz-1; k++)
      {
        for(int i4=1; i4<nt-1; i4++)
        {
          row = S(i,j,k,i4,ny,nz,nt); // equation 
          // ccxx*uxx+ccyy*uyy+cczz*uzz+cctt*utt+
          // ccxy*uxy+ccxz*uxz+ccxt*uxt+ccyz*uyz++ccyt*uyt++cczt*uzt+
          // ccx*ux+ccy*uy+ccz*uz+cct*ut+ccc*u=f(x,y,z)
          // for(each term 

          // ccxx*uxx  =  d^2 uana(x,y,z,t)/dx^2 
          new nuderiv(2,nx,i,xg,cxx);
          for(int ki=0; ki<nx; ki++)
          {
	    // if(i==1 && j==1) System.out.print("cxx["+ki+"]="+cxx[ki]+" \n");
            if(ki==0 || ki==nx-1)
            {
              ut[row][cs] = ut[row][cs] - ccxx*cxx[ki]*
                            uana(xg[ki],yg[j],zg[k],tg[i4]);
            }	 
            else
            {
              ut[row][S(ki,j,k,i4,ny,nz,nt)] = ut[row][S(ki,j,k,i4,ny,nz,nt)] +
                                               ccxx*cxx[ki];
            }
          } // ki

          // ccyy*uyy  = d^2 uana(x,y,z,t)/dy^2 
          new nuderiv(2,ny,j,yg,cyy);
          for(int kj=0; kj<ny; kj++)
          {
	    // if(i==1 && j==1) System.out.print("cyy["+kj+"]="+cyy[kj]+" \n");
            if(kj==0 || kj==ny-1)
	    {
	      ut[row][cs] = ut[row][cs] - ccyy*cyy[kj]*
                            uana(xg[i],yg[kj],zg[k],tg[i4]);
            }
            else
	    {
	      ut[row][S(i,kj,k,i4,ny,nz,nt)] = ut[row][S(i,kj,k,i4,ny,nz,nt)] +
                                               ccyy*cyy[kj];
	    }
          } // kj

          // cczz*uzz  = d^2 uana(x,y,z,t)/dz^2 
          new nuderiv(2,nz,k,zg,czz);
          for(int kk=0; kk<nz; kk++)
          {
	    // if(i==1 && j==1) System.out.print("czz["+kk+"]="+czz[kk]+" \n");
            if(kk==0 || kk==nz-1)
	    {
	      ut[row][cs] = ut[row][cs] - cczz*czz[kk]*
                            uana(xg[i],yg[j],zg[kk],tg[i4]);
            }
            else
	    {
	      ut[row][S(i,j,kk,i4,ny,nz,nt)] = ut[row][S(i,j,kk,i4,ny,nz,nt)] +
                                               cczz*czz[kk];
	    }
          } // kk

          // cctt*utt  = d^2 uana(x,y,z,t)/dt^2 
          new nuderiv(2,nt,i4,tg,ctt);
          for(int j4=0; j4<nt; j4++)
          {
	    // if(i==1 && j==1) System.out.print("ctt["+j4+"]="+ctt[j4]+" \n");
            if(j4==0 || j4==nt-1)
	    {
	      ut[row][cs] = ut[row][cs] - cctt*ctt[j4]*
                            uana(xg[i],yg[j],zg[k],tg[j4]);
            }
            else
	    {
	      ut[row][S(i,j,k,j4,ny,nz,nt)] = ut[row][S(i,j,k,j4,ny,nz,nt)] +
                                               cctt*ctt[j4];
	    }
          } // j4

          // ccxy*uxy  =  d^2 uana(x,y,z,t)/dxdy 
          new nuderiv(1,nx,i,xg,cx);
          new nuderiv(1,ny,j,yg,cy);
          for(int ki=0; ki<nx; ki++)
          {
            for(int kj=0; kj<ny; kj++)
            {
              if(ki==0 || ki==nx-1 || kj==0 || kj==ny-1)
              {
	        ut[row][cs] = ut[row][cs] - ccxy*cx[ki]*cy[kj]*
                              uana(xg[ki],yg[kj],zg[k],tg[i4]);
              }
              else
	      {
	        ut[row][S(ki,kj,k,i4,ny,nz,nt)] = ut[row][S(ki,kj,k,i4,ny,nz,nt)] +
                                                  ccxy*cx[ki]*cy[kj];
	      }
            } // kj 
          } // ki 

          // ccxz*uxz  =  d^2 uana(x,y,z,t)/dxdz 
          new nuderiv(1,nx,i,xg,cx);
          new nuderiv(1,nz,k,zg,cz);
          for(int ki=0; ki<nx; ki++)
          {
            for(int kk=0; kk<nz; kk++)
            {
              if(ki==0 || ki==nx-1 || kk==0 || kk==nz-1)
              {
	        ut[row][cs] = ut[row][cs] - ccxz*cx[ki]*cz[kk]*
                                            uana(xg[ki],yg[j],zg[kk],tg[i4]);
              }
              else
	      {
	        ut[row][S(ki,j,kk,i4,ny,nz,nt)] = ut[row][S(ki,j,kk,i4,ny,nz,nt)] +
                                                  ccxz*cx[ki]*cz[kk];
	      }
            } // kk 
          } // ki 

          // ccxt*uxt  =  d^2 uana(x,y,z,t)/dxdt 
          new nuderiv(1,nx,i,xg,cx);
          new nuderiv(1,nt,i4,tg,ct);
          for(int ki=0; ki<nx; ki++)
          {
            for(int j4=0; j4<nt; j4++)
            {
              if(ki==0 || ki==nx-1 || j4==0 || j4==nt-1)
              {
	        ut[row][cs] = ut[row][cs] - ccxt*cx[ki]*ct[j4]*
                                            uana(xg[ki],yg[j],zg[k],tg[j4]);
              }
              else
	      {
	        ut[row][S(ki,j,k,j4,ny,nz,nt)] = ut[row][S(ki,j,k,j4,ny,nz,nt)] +
                                                  ccxt*cx[ki]*ct[j4];
	      }
            } // j4 
          } // ki

          // ccyz*uyz  =  d^2 uana(x,y,z,t)/dydz 
          new nuderiv(1,ny,j,yg,cy);
          new nuderiv(1,nz,k,zg,cz);
          for(int kj=0; kj<ny; kj++)
          {
            for(int kk=0; kk<nz; kk++)
            {
              if(kj==0 || kj==ny-1 || kk==0 || kk==nz-1)
              {
	        ut[row][cs] = ut[row][cs] - ccyz*cy[kj]*cz[kk]*
                              uana(xg[i],yg[kj],zg[kk],tg[i4]);
              }
              else
	      {
	        ut[row][S(i,kj,kk,i4,ny,nz,nt)] = ut[row][S(i,kj,kk,i4,ny,nz,nt)] +
                                                  ccyz*cy[kj]*cz[kk];
	      }
            } // kk
          } // kj 

          // ccyt*uyt  =  d^2 uana(x,y,z,t)/dydt 
          new nuderiv(1,ny,j,yg,cy);
          new nuderiv(1,nt,i4,tg,ct);
          for(int kj=0; kj<ny; kj++)
          {
            for(int j4=0; j4<nt; j4++)
            {
              if(kj==0 || kj==ny-1 || j4==0 || j4==nt-1)
              {
	        ut[row][cs] = ut[row][cs] - ccyt*cy[kj]*ct[j4]*
                              uana(xg[i],yg[kj],zg[k],tg[j4]);
              }
              else
	      {
	        ut[row][S(i,kj,k,j4,ny,nz,nt)] = ut[row][S(i,kj,k,j4,ny,nz,nt)] +
                                                  ccyt*cy[kj]*ct[j4];
	      }
            } // j4
          } // kj 

          // cczt*uzt  =  d^2 uana(x,y,z,t)/dzdt 
          new nuderiv(1,nz,k,zg,cz);
          new nuderiv(1,nt,i4,tg,ct);
          for(int kk=0; kk<nz; kk++)
          {
            for(int j4=0; j4<nt; j4++)
            {
              if(kk==0 || kk==nz-1 || j4==0 || j4==nt-1)
              {
	        ut[row][cs] = ut[row][cs] - cczt*cz[kk]*ct[j4]*
                              uana(xg[i],yg[j],zg[kk],tg[j4]);
              }
              else
	      {
	        ut[row][S(i,j,kk,j4,ny,nz,nt)] = ut[row][S(i,j,kk,j4,ny,nz,nt)] +
                                                  cczt*cz[kk]*ct[j4];
	      }
            } // kk
          } // j4


          // ccx*ux  =  d uana(x,y,z,t)/dx 
          new nuderiv(1,nx,i,xg,cx);
          for(int ki=0; ki<nx; ki++)
          {
	    // if(i==1 && j==1) System.out.print("cx["+ki+"]="+cx[ki]+" \n");
            if(ki==0 || ki==nx-1)
            {
	      ut[row][cs] = ut[row][cs] - ccx*cx[ki]*
                            uana(xg[ki],yg[j],zg[k],tg[i4]);
            }	 
            else
	    {
	      ut[row][S(ki,j,k,i4,ny,nz,nt)] = ut[row][S(ki,j,k,i4,ny,nz,nt)] +
                                               ccx*cx[ki];
	    }
          } // ki

          // ccy*uy  = d uana(x,y,z,t)/dy 
          new nuderiv(1,ny,j,yg,cy);
          for(int kj=0; kj<ny; kj++)
          {
	    // if(i==1 && j==1) System.out.print("cy["+kj+"]="+cy[kj]+" \n");
            if(kj==0 || kj==ny-1)
            {
	      ut[row][cs] = ut[row][cs] - ccy*cy[kj]*
                            uana(xg[i],yg[kj],zg[k],tg[i4]);
            }
            else
            {
	      ut[row][S(i,kj,k,i4,ny,nz,nt)] = ut[row][S(i,kj,k,i4,ny,nz,nt)] +
                                               ccy*cy[kj];
            }
          } // kj

          // ccz*uz  = d uana(x,y,z,t)/dz 
          new nuderiv(1,nz,k,zg,cz);
          for(int kk=0; kk<nz; kk++)
          {
	    // if(i==1 && j==1) System.out.print("cz["+kk+"]="+cz[kk]+" \n");
            if(kk==0 || kk==nz-1)
            {
	      ut[row][cs] = ut[row][cs] - ccz*cz[kk]*
                            uana(xg[i],yg[j],zg[kk],tg[i4]);
            }
            else
            {
	      ut[row][S(i,j,kk,i4,ny,nz,nt)] = ut[row][S(i,j,kk,i4,ny,nz,nt)] +
                                               ccz*cz[kk];
            }
          } // kk

          // cct*ut  = d uana(x,y,z,t)/dt 
          new nuderiv(1,nt,i4,tg,ct);
          for(int j4=0; j4<nt; j4++)
          {
	    // if(i==1 && j==1) System.out.print("ct["+j4+"]="+ct[j4]+" \n");
            if(j4==0 || j4==nt-1)
            {
	      ut[row][cs] = ut[row][cs] - cct*ct[j4]*
                            uana(xg[i],yg[j],zg[k],tg[j4]);
            }
            else
            {
	      ut[row][S(i,j,k,j4,ny,nz,nt)] = ut[row][S(i,j,k,j4,ny,nz,nt)] +
                                               cct*ct[j4];
            }
          } // j4


          // ccc*u 
          ut[row][S(i,j,k,i4,ny,nz,nt)] = ut[row][S(i,j,k,i4,ny,nz,nt)] + ccc;


          // RHS constant 
          ut[row][cs] = ut[row][cs] + f(xg[i],yg[j],zg[k],tg[i4]);
          // end terms
        } // i4 
      } // k
    } // j 
  } // i 
} // end init_matrix 

void check_soln(int nx, int ny, int nz, int nt, double us[],
                double xg[], double yg[], double zg[], double tg[],
                double ccxx, double ccyy, double cczz, double cctt,
                double ccxy, double ccxz, double ccxt,
                double ccyz, double ccyt, double cczt,
                double ccx, double ccy, double ccz, double cct, double ccc)
                // check when solution unknown
{
  double cxx[] = new double[10];  // new nuderiv coefficients 
  double cyy[] = new double[10];
  double czz[] = new double[10];
  double ctt[] = new double[10];
  double cx[] = new double[10]; 
  double cy[] = new double[10];
  double cz[] = new double[10];
  double ct[] = new double[10];
  int row;
  int cs = (nx-2)*(ny-2)*(nz-2)*(nt-2);
  double smaxerr, err, uxx=0.0, uyy=0.0, uzz=0.0, utt=0.0;
  double uxy=0.0, uxz=0.0, uxt=0.0, uyz=0.0, uyt=0.0, uzt=0.0;
  double ux=0.0, uy=0.0, uz=0.0, ut=0.0;

  // check ccxx*uxx(x,y,z,t)+ccyy*uyy(x,y,z,t)+cczz*uzz(x,y,z,t)+
  //       cctt*utt(x,y,z,t)+ccxy*uxy(x,y,z,t)+ccxy*uxy(x,y,z,t)+
  //       ccxt*uxt(x,y,z,t)+ccyz*uyz(x,y,z,t)+ccyt*uyz(x,y,z,t)+
  //       cczt*uzt(x,y,z,t)+ccx * ux(x,y,z,t)+ccy * uy(x,y,z,t)+
  //       ccz * uz(x,y,z,t)+cct * ut(x,y,z,t)+ccc*uana(x,y,z,t)=f(x,y,z,t)
  // against computed solution us[S(i,j,k,i4,ny,nz,nt)] at each i,j,k,i4 
  System.out.print("check computed solution against PDE \n");
  smaxerr = 0.0;

  for(int i=1; i<nx-1; i++)  // through dof equations 
  {
    new nuderiv(1,nx,i,xg,cx); // get discrete derivative coefficients
    new nuderiv(2,nx,i,xg,cxx);
    for(int j=1; j<ny-1; j++)
    {
      new nuderiv(1,ny,j,yg,cy);
      new nuderiv(2,ny,j,yg,cyy);
      for(int k=1; k<nz-1; k++)
      {
        new nuderiv(1,nz,k,zg,cz);
        new nuderiv(2,nz,k,zg,czz);
        for(int i4=1; i4<nt-1; i4++)
        {
          new nuderiv(1,nt,i4,tg,ct);
          new nuderiv(2,nt,i4,tg,ctt);

          ux  = 0.0;
          uxx = 0.0;
          for(int ki=0; ki<nx; ki++) // compute numerical derivative 
          {
            if(ki==0 || ki==nx-1)
            {
              ux  = ux  + cx [ki]*uana(xg[ki],yg[j],zg[k],tg[i4]); // boundary 
              uxx = uxx + cxx[ki]*uana(xg[ki],yg[j],zg[k],tg[i4]); // boundary 
            }
            else
            {
              ux  = ux  + cx [ki]*us[S(ki,j,k,i4,ny,nz,nt)];
              uxx = uxx + cxx[ki]*us[S(ki,j,k,i4,ny,nz,nt)];
            }
          } // ki

          uy  = 0.0;
          uyy = 0.0;
          for(int kj=0; kj<ny; kj++)
          {
            if(kj==0 || kj==ny-1)
            {
              uy  = uy  + cy [kj]*uana(xg[i],yg[kj],zg[k],tg[i4]); // boundary 
              uyy = uyy + cyy[kj]*uana(xg[i],yg[kj],zg[k],tg[i4]); // boundary 
            }
            else
            {
              uy  = uy  + cy [kj]*us[S(i,kj,k,i4,ny,nz,nt)];
              uyy = uyy + cyy[kj]*us[S(i,kj,k,i4,ny,nz,nt)];
            }
          } // kj

          uz  = 0.0;
          uzz = 0.0;
          for(int kk=0; kk<nz; kk++)
          {
            if(kk==0 || kk==nz-1)
            {
              uz  = uz  + cz [kk]*uana(xg[i],yg[j],zg[kk],tg[i4]); // boundary 
              uzz = uzz + czz[kk]*uana(xg[i],yg[j],zg[kk],tg[i4]); // boundary 
            }
            else
            {
              uz  = uz  + cz [kk]*us[S(i,j,kk,i4,ny,nz,nt)];
              uzz = uzz + czz[kk]*us[S(i,j,kk,i4,ny,nz,nt)];
            }
          } // kk

          ut  = 0.0;
          utt = 0.0;
          for(int j4=0; j4<nt; j4++)
          {
            if(j4==0 || j4==nt-1)
            {
              ut  = ut  + ct [j4]*uana(xg[i],yg[j],zg[k],tg[j4]); // boundary 
              utt = utt + ctt[j4]*uana(xg[i],yg[j],zg[k],tg[j4]); // boundary 
            }
            else
            {
              ut  = ut  + ct [j4]*us[S(i,j,k,j4,ny,nz,nt)];
              utt = utt + ctt[j4]*us[S(i,j,k,j4,ny,nz,nt)];
            }
          } // j4

          uxy = 0.0;
          for(int ki=0; ki<nx; ki++)
          {
            for(int kj=0; kj<ny; kj++)
            {
              if(ki==0 || ki==nx-1 || kj==0 || kj==ny-1)
              {
                uxy  = uxy  + cx[ki]*cy[kj]*uana(xg[ki],yg[kj],zg[k],tg[i4]); // bound
              }
              else
              {
                uxy = uxy + cx[ki]*cy[kj]*us[S(ki,kj,k,i4,ny,nz,nt)];
              }
            } // kj
          } // ki

          uxz = 0.0;
          for(int ki=0; ki<nx; ki++)
          {
            for(int kk=0; kk<nz; kk++)
            {
              if(ki==0 || ki==nx-1 || kk==0 || kk==nz-1)
              {
                uxz  = uxz  + cx[ki]*cz[kk]*uana(xg[ki],yg[j],zg[kk],tg[i4]); // bound
              }
              else
              {
                uxz = uxz + cx[ki]*cz[kk]*us[S(ki,j,kk,i4,ny,nz,nt)];
              }
            } // kk
          } // ki

          uxt = 0.0;
          for(int ki=0; ki<nx; ki++)
          {
            for(int j4=0; j4<nt; j4++)
            {
              if(ki==0 || ki==nx-1 || j4==0 || j4==nt-1)
              {
                uxt  = uxt  + cx[ki]*ct[j4]*uana(xg[ki],yg[j],zg[k],tg[j4]); // bound
              }
              else
              {
                uxt = uxt + cx[ki]*ct[j4]*us[S(ki,j,k,j4,ny,nz,nt)];
              }
            } // j4
          } // ki

          uyz = 0.0;
          for(int kj=0; kj<ny; kj++)
          {
            for(int kk=0; kk<nz; kk++)
            {
              if(kj==0 || kj==ny-1 || kk==0 || kk==nz-1)
              {
                uyz  = uyz  + cy[kj]*cz[kk]*uana(xg[i],yg[kj],zg[kk],tg[i4]); // bound
              }
              else
              {
                uyz = uyz + cy[kj]*cz[kk]*us[S(i,kj,kk,i4,ny,nz,nt)];
              }
            } // kk
          } // kj

          uyt = 0.0;
          for(int kj=0; kj<ny; kj++)
          {
            for(int j4=0; j4<nt; j4++)
            {
              if(kj==0 || kj==ny-1 || j4==0 || j4==nt-1)
              {
                uyt  = uyt  + cy[kj]*ct[j4]*uana(xg[i],yg[kj],zg[k],tg[j4]); // bound
              }
              else
              {
                uyt = uyt + cy[kj]*ct[j4]*us[S(i,kj,k,j4,ny,nz,nt)];
              }
            } // j4
          } // kj

          uzt = 0.0;
          for(int kk=0; kk<nz; kk++)
          {
            for(int j4=0; j4<nt; j4++)
            {
              if(kk==0 || kk==nz-1 || j4==0 || j4==nt-1)
              {
                uzt  = uzt  + cz[kk]*ct[j4]*uana(xg[i],yg[j],zg[kk],tg[j4]); // bound
              }
              else
              {
                uzt = uzt + cz[kk]*ct[j4]*us[S(i,j,kk,j4,ny,nz,nt)];
              }
            } // j4
          } // kk

          err = ccxx*uxx + ccyy*uyy + cczz*uzz + cctt*utt +
                ccxy*uxy + ccxz*uxz + ccxt*uxt + ccyz*uyz +
                ccyt*uyt + cczt*uzt + ccx*ux   + ccy*uy   + ccz*uz   +
                cct*ut   + ccc*us[S(i,j,k,i4,ny,nz,nt)]
              - f(xg[i],yg[j],zg[k],tg[i4]); // minus right hand side 
          if(Math.abs(err) > smaxerr) smaxerr = Math.abs(err);
        } // i4
      } // k z
    } // j y 
  } // i x 

  System.out.print(" \n");
  System.out.print("check computed soln against PDE max error="+smaxerr+" \n");
  System.out.print(" \n");
  System.out.print("x="+xg[nx-2]+", y="+yg[ny-2]+", z="+zg[ny-2]+" \n");
  System.out.print("t="+tg[nt-2]+", uxx="+uxx+", uyy="+uyy+", uzz="+uzz+" \n");
  System.out.print("utt="+utt+", ux="+ux+", uy="+uy+", uz="+uz+", ut="+ut+" \n");
  System.out.print("us="+us[S(nx-2,ny-2,nz-2,nt-2,ny,nz,nt)]+" \n");
  System.out.print("f(x,y,z,t)="+f(xg[nx-2],yg[ny-2],zg[nz-2],tg[nt-2])+" \n");
  System.out.print(" \n");
} // end Check_Soln 

double uana(double x,double y,double z,double t) // ub(x,y,z,t) for boundary 
{
  return x*x+y*y+z*z+t*t+x*y+x*z+x*t+y*z+y*t+z*t+1.0;
} // end uana

double f(double x,double y,double z,double t) // RHS f(x,y,z,t)
{
    return 15.0*(x*x+y*y+z*z+t*t+x*y+x*z+x*t+y*z+y*t+z*t)+
           61.0*x+62.0*y+63.0*z+64.0*t+80.0;
} // end f

  public static void main (String[] args)
  {
    new pde24js_eq();
  }
} // end pde24js_eq.java