// pde46h_eq.java  homogeneous Biharmonic PDE in six dimensions
// solve Uvvvv(v,w,x,y,z,t) + Uwwww(v,w,x,y,z,t) + Uxxxx(v,w,x,y,z,t) + 
//       Uyyyy(v,w,x,y,z,t) + Uzzzz(v,w,x,y,z,t) + Utttt(v,w,x,y,z,t) + 
//       2 Uvv(v,w,x,y,z,t) + 2 Uww(v,w,x,y,z,t) + 2 Uxx(v,w,x,y,z,t) +
//       2 Uyy(v,w,x,y,z,t) + 2 Uzz(v,w,x,y,z,t) + 2 Utt(v,w,x,y,z,t) + 
//         6 U(v,w,x,y,z,t) = 0
//
// boundary conditions computed using uana(v,w,x,y,z,t)
// analytic solution is U(v,w,x,y,z,t) = sin(v+w+x+y+z+t)
//
// subscripts in code use v[i], w[ii], x[iii], y[iiii], z[iiiii], t[iiiiii]
//
// replace continuous derivatives with discrete derivatives
// solve for U(w,x,y,z,t) numerical approximation U
// A * U = F, solve simultaneous equations for U
//
 
import java.text.*;
import java.io.*;
// uses rderiv.java
// uses simeq.java

class pde46h_eq
{
  int nv = 6;
  int nw = 6;
  int nx = 6;
  int ny = 6;
  int nz = 6;
  int nt = 6;
  int nvwxyzt = (nv-2)*(nw-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2); // DOF
  double A[] = new double[nvwxyzt*(nvwxyzt+1)]; // last column is RHS 
  double U[] = new double[nvwxyzt]; // comupted solution, no boundary
  double vg[] = new double[nv]; // coordinates
  double wg[] = new double[nw];
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double tg[] = new double[nt];
  double Ua[] = new double[nvwxyzt]; // analytic solution for checking, no bou
  double vmin = 0.0; // problem parameters 
  double vmax = 1.5708;
  double wmin = 0.0;
  double wmax = 1.5708;
  double xmin = 0.0;
  double xmax = 1.5708;
  double ymin = 0.0;
  double ymax = 1.5708;
  double zmin = 0.0;
  double zmax = 1.5708;
  double tmin = 0.0;
  double tmax = 1.5708;
  double hv, hw, hx, hy, hz, ht;

  pde46h_eq()
  {
    double v, w, x, y, z, t;
    int i,ii,iii,iiii,iiiii,iiiiii,j,jj,jjj,jjjj,jjjjj,jjjjjj;
    int k, cs;
    // coeficients of terms in first equation 
    // derivative point arrays, dynamic in this version 
    double cvv[] = new double[nv];
    double cww[] = new double[nw];
    double cxx[] = new double[nx];
    double cyy[] = new double[ny];
    double czz[] = new double[nz];
    double ctt[] = new double[nt]; 
    double cvvvv[] = new double[nv];
    double cwwww[] = new double[nw];
    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("pde46h_eq.java running ");
    System.out.println("Given uvvvv(v,w,x,y,z,t) + uwwww(v,w,x,y,z,t) + uxxxx(v,w,x,y,z,t) +");
    System.out.println("      uyyyy(v,w,x,y,z,t) + uzzzz(v,w,x,y,z,t) + uzzzz(v,w,x,y,z,t) +");
    System.out.println("      2 uvv(v,w,x,y,z,t) + 2 uww(v,w,x,y,z,t) + 2 uxx(v,w,x,y,z,t) +");
    System.out.println("      2 uyy(v,w,x,y,z,t) + 2 uzz(v,w,x,y,z,t) + 2 utt(v,w,x,y,z,t) +");
    System.out.println("      6   u(v,w,x,y,z,t) = 0 ");
    System.out.println(" "); 
    System.out.println("vmin<=v<=vmax wmin<=w<=wmax xmin<=x<=xmax ");
    System.out.println("ymin<=y<=ymax zmin<=z<=zmax tmin<=t<=tmax Boundaries ");
    System.out.println("Analytic solution U(v,w,x,y,z,t) =  sin(v+w+x+y+z+t) ");
    System.out.println(" ");

    System.out.println("vmin="+vmin+", vmax="+vmax+
                       ", wmin="+wmin+", wmax="+wmax);
    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("nv="+nv+", nw="+nw+", nx="+nx+", ny="+ny+
                     ", nz="+nz+", nt="+nt);
  time_start = (double)System.currentTimeMillis()/1000.0;
  System.out.println("time start="+time_start+" seconds");
  System.out.println("DOF, degrees of freedom, nvwxyzt="+nvwxyzt);

  time_last = time_start;

  System.out.println(" ");
  System.out.println("v grid and analytic solution at minimums ");
  hv = (vmax-vmin)/(double)(nv-1);
  for(i=0; i<nv; i++)
  {
    vg[i] = vmin+(double)i*hv;
    System.out.println("i="+i+", Ua("+vg[i]+","+wmin+","+xmin+
                       ","+ymin+","+zmin+","+tmin+")="+
                       uana(vg[i],wmin,xmin,ymin,zmin,tmin)); 
  } // i   
  System.out.println(" ");
  System.out.println("w grid and analytic solution at minimums ");
  hw = (wmax-wmin)/(double)(nw-1);
  for(ii=0; ii<nw; ii++)
  {
    wg[ii] = wmin+(double)ii*hw;
    System.out.println("ii="+ii+", Ua("+vmin+","+wg[ii]+","+xmin+
                       ","+ymin+","+zmin+","+tmin+")="+
                       uana(vmin,wg[ii],xmin,ymin,zmin,tmin)); 
  } // ii   
  System.out.println(" ");
  System.out.println("x grid and analytic solution at minimums ");
  hx = (xmax-xmin)/(double)(nx-1);
  for(iii=0; iii<nx; iii++)
  {
    xg[iii] = xmin+(double)iii*hx;
    System.out.println("iii="+iii+", Ua("+vmin+","+wmin+","+xg[iii]+
                       ","+ymin+","+zmin+","+tmin+")="+
                       uana(vmin,wmin,xg[iii],ymin,zmin,tmin)); 
  } // iii   
  System.out.println(" ");
  System.out.println("y grid and analytic solution at minimums ");
  hy = (ymax-ymin)/(double)(ny-1);
  for(iiii=0; iiii<ny; iiii++)
  {
    yg[iiii] = ymin+(double)iiii*hy;
    System.out.println("iiii="+iiii+", Ua("+vmin+","+wmin+","+xmin+
                       ","+yg[iiii]+","+zmin+","+tmin+")="+
                       uana(vmin,wmin,xmin,yg[iiii],zmin,tmin)); 
  } // iiii 
  System.out.println(" ");
  System.out.println("z grid and analytic solution at minimums ");
  hz = (zmax-zmin)/(double)(nz-1);
  for(iiiii=0; iiiii<nz; iiiii++)
  {
    zg[iiiii] = zmin+(double)iiiii*hz;
    System.out.println("iiiii="+iiiii+", Ua("+vmin+","+wmin+","+xmin+
                       ","+ymin+","+zg[iiiii]+","+tmin+")="+
                       uana(vmin,wmin,xmin,ymin,zg[iiiii],tmin)); 
  } // iiiii 
  System.out.println(" ");
  System.out.println("t grid and analytic solution at minimiums ");
  ht = (tmax-tmin)/(double)(nt-1);
  for(iiiiii=0; iiiiii<nt; iiiiii++)
  {
    tg[iiiiii] = tmin+(double)iiiiii*ht;
    System.out.println("iiiiii="+iiiiii+", Ua("+vmin+","+wmin+","+xmin+
                       ","+ymin+","+zmin+","+tg[iiiiii]+")="+
                       uana(vmin,wmin,xmin,ymin,zmin,tg[iiiiii])); 
  } // iiiiii 
  System.out.println(" ");

  System.out.println("check rderiv(4,nv,0,hv,cvvvv ");
  new rderiv(4,nv,0,hv,cvvvv);
  for(i=0; i<nv; i++) System.out.println("cvvvv["+i+"]="+cvvvv[i]);
  System.out.println("check rderiv(4,nv,nv-1,hv,cvvvv ");
  new rderiv(4,nv,nv-1,hv,cvvvv);
  for(i=0; i<nv; i++) System.out.println("cvvvv["+i+"]="+cvvvv[i]);
  System.out.println("");

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

  System.out.println("initialize A ");
  // initialize A 
  for(i=1; i<nv-1; i++)
  {
    for(ii=1; ii<nw-1; ii++)
    {
      for(iii=1; iii<nx-1; iii++)
      {
        for(iiii=1; iiii<ny-1; iiii++)
        {
          for(iiiii=1; iiiii<nz-1; iiiii++)
          {
            for(iiiiii=1; iiiiii<nt-1; iiiiii++)
            {
              for(j=1; j<nv-1; j++)
              {
                for(jj=1; jj<nw-1; jj++)
                {
                  for(jjj=1; jjj<nx-1; jjj++)
                  {
                    for(jjjj=1; jjjj<ny-1; jjjj++)
                    {
                      for(jjjjj=1; jjjjj<nz-1; jjjjj++)
                      {
                        for(jjjjjj=1; jjjjjj<nt-1; jjjjjj++)
                        {
                          A[ss(i,ii,iii,iiii,iiiii,iiiiii,
                               j,jj,jjj,jjjj,jjjjj,jjjjjj)] = 0.0;
	                } // jjjjjj 
	              } // jjjjj 
	            } // jjjj 
	          } // jjj 
	        } // jj 
              } // j 
            } // iiiiii 
          } // iiiii 
        } // iiii 
      } // iii 
    } // ii 
  } // i 

  time_now = (double)System.currentTimeMillis()/1000.0;
  System.out.println("time now ="+(time_now-time_start)+
                     " seconds, 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 = nvwxyzt; // constant RHS column 
  val = 0.0;
  for(i=1; i<nv-1; i++)
    {
    for(ii=1; ii<nw-1; ii++)
    {
      for(iii=1; iii<nx-1; iii++)
      {
        for(iiii=1; iiii<ny-1; iiii++)
        {
          for(iiiii=1; iiiii<nz-1; iiiii++)
          {
            for(iiiiii=1; iiiiii<nt-1; iiiiii++)
            { 
              k = s(i,ii,iii,iiii,iiiii,iiiiii); // row index 

              // for each term in first equation 

	      // for d^4u/dv^4 
              new rderiv(4,nv,i,hv,cvvvv);
              for(j=0; j<nv; j++)
              {
                val = cvvvv[j];
                if(j==0 || j==nv-1)
	          A[k*(nvwxyzt+1)+cs] -= val*
		    uana(vg[j],wg[ii],xg[iii],yg[iiii],zg[iiiii],tg[iiiiii]);
                else A[sk(k,j,ii,iii,iiii,iiiii,iiiiii)] += val;
	      }

	      // for d^4u/dw^4 
              new rderiv(4,nw,ii,hw,cwwww);
              for(jj=0; jj<nw; jj++)
              {
                val = cwwww[jj];
                if(jj==0 || jj==nw-1)
	          A[k*(nvwxyzt+1)+cs] -= val*
		    uana(vg[i],wg[jj],xg[iii],yg[iiii],zg[iiiii],tg[iiiiii]);
                else A[sk(k,i,jj,iii,iiii,iiiii,iiiiii)] += val;
	      }

	      // for d^4u/dx^4 
              new rderiv(4,nx,iii,hx,cxxxx);
              for(jjj=0; jjj<nx; jjj++)
              {
                val = cxxxx[jjj];
                if(jjj==0 || jjj==nx-1)
		  A[k*(nvwxyzt+1)+cs] -= val*
		    uana(vg[i],wg[ii],xg[jjj],yg[iiii],zg[iiiii],tg[iiiiii]);
                else A[sk(k,i,ii,jjj,iiii,iiiii,iiiiii)] += val;
	      }

	      // for d^4u/dy^4 
              new rderiv(4,ny,iiii,hy,cyyyy);
              for(jjjj=0; jjjj<ny; jjjj++)
              {
                val = cyyyy[jjjj];
                if(jjjj==0 || jjjj==ny-1)
	          A[k*(nvwxyzt+1)+cs] -=  val*
		    uana(vg[i],wg[ii],xg[iii],yg[jjjj],zg[iiiii],tg[iiiiii]);
                else A[sk(k,i,ii,iii,jjjj,iiiii,iiiiii)] += val;
	      }

	      // for d^4u/dz^4 
              new rderiv(4,nz,iiiii,hz,czzzz);
              for(jjjjj=0; jjjjj<nz; jjjjj++)
              {
                val = czzzz[jjjjj];
                if(jjjjj==0 || jjjjj==nz-1)
	          A[k*(nvwxyzt+1)+cs] -= val*
		    uana(vg[i],wg[ii],xg[iii],yg[iiii],zg[jjjjj],tg[iiiiii]);
                else A[sk(k,i,ii,iii,iiii,jjjjj,iiiiii)] += val;
	      }

	      // for d^4u/dt^4 
              new rderiv(4,nt,iiiiii,ht,ctttt);
              for(jjjjjj=0; jjjjjj<nt; jjjjjj++)
              {
                val = ctttt[jjjjjj];
                if(jjjjjj==0 || jjjjjj==nt-1)
	          A[k*(nvwxyzt+1)+cs] -= val*
		    uana(vg[i],wg[ii],xg[iii],yg[iiii],zg[iiiii],tg[jjjjjj]);
                else A[sk(k,i,ii,iii,iiii,iiiii,jjjjjj)] += val;
	      }


	      // for 2 d^2u/dv^2 
              new rderiv(2,nv,ii,hv,cvv);
              for(j=0; j<nv; j++)
              {
                val = 2.0*cvv[j];
                if(j==0 || j==nv-1)
	          A[k*(nvwxyzt+1)+cs] -= val*
		    uana(vg[j],wg[ii],xg[iii],yg[iiii],zg[iiiii],tg[iiiiii]);
                else A[sk(k,j,ii,iii,iiii,iiiii,iiiiii)] += val;
	      }

	      // for 2 d^2u/dw^2 
              new rderiv(2,nw,ii,hw,cww);
              for(jj=0; jj<nw; jj++)
              {
                val = 2.0*cww[jj];
                if(jj==0 || jj==nw-1)
	          A[k*(nvwxyzt+1)+cs] -= val*
		    uana(vg[i],wg[jj],xg[iii],yg[iiii],zg[iiiii],tg[iiiiii]);
                else A[sk(k,i,jj,iii,iiii,iiiii,iiiiii)] += val;
	      }

	      // for 2 d^2u/dx^2 
              new rderiv(2,nx,iii,hx,cxx);
              for(jjj=0; jjj<nx; jjj++)
              {
                val = 2.0*cxx[jjj];
                if(jjj==0 || jjj==nx-1)
	          A[k*(nvwxyzt+1)+cs] -= val*
		    uana(vg[i],wg[ii],xg[jjj],yg[iiii],zg[iiiii],tg[iiiiii]);
                else A[sk(k,i,ii,jjj,iiii,iiiii,iiiiii)] += val;
	      }

	      // for 2 d^2u/dy^2 
              new rderiv(2,ny,iiii,hy,cyy);
              for(jjjj=0; jjjj<ny; jjjj++)
              {
                val = 2.0*cyy[jjjj];
                if(jjjj==0 || jjjj==ny-1)
	          A[k*(nvwxyzt+1)+cs] -=  val*
		    uana(vg[i],wg[ii],xg[iii],yg[jjjj],zg[iiiii],tg[iiiiii]);
                else A[sk(k,i,ii,iii,jjjj,iiiii,iiiiii)] += val;
	      }

	      // for 2 d^2u/dz^2 
              new rderiv(2,nz,iiiii,hz,czz);
              for(jjjjj=0; jjjjj<nz; jjjjj++)
              { 
                val = 2.0*czz[jjjjj];
                if(jjjjj==0 || jjjjj==nz-1)
	          A[k*(nvwxyzt+1)+cs] -= val*
		    uana(vg[i],wg[ii],xg[iii],yg[iiii],zg[jjjjj],tg[iiiiii]);
                else A[sk(k,i,ii,iii,iiii,jjjjj,iiiiii)] += val;
	      }

	      // for 2 d^2u/dt^2 
              new rderiv(2,nt,iiiiii,ht,ctt);
              for(jjjjjj=0; jjjjjj<nt; jjjjjj++)
              {
                val = 2.0*ctt[jjjjjj];
                if(jjjjjj==0 || jjjjjj==nt-1)
	          A[k*(nvwxyzt+1)+cs] -= val*
		    uana(vg[i],wg[ii],xg[iii],yg[iiii],zg[iiiii],tg[jjjjjj]);
                else A[sk(k,i,ii,iii,iiii,iiiii,jjjjjj)] += val;
	      }

	      // for 6 u(v,w,x,y,z,t) 
              A[sk(k,i,ii,iii,iiii,iiiii,iiiiii)] += 6.0;

              // for f(v,w,x,y,z,t) RHS zero 

              // end terms 

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

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

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

    new simeq(nvwxyzt, A, U); // A destroyed 

    System.out.println("equations solved ");
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("time now = "+(time_now-time_start)+
                       " seconds, time for previous section ="+
                       (time_now-time_last)+" seconds");
    time_last = time_now;
    System.out.println("");
  
    System.out.println("check PDE against computed solution");
    check_soln(U);
    write_soln(U);
    System.out.println("");

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("          U computed, Ua analytic, error ");
    for(i=1; i<nv-1; i++)
    {
      for(ii=1; ii<nw-1; ii++)
      {
        for(iii=1; iii<nx-1; iii++)
        {
          for(iiii=1; iiii<ny-1; iiii++)
          {
            for(iiiii=1; iiiii<nz-1; iiiii++)
            {
              for(iiiiii=1; iiiiii<nt-1; iiiiii++)
              {
                err = Math.abs(U[s(i,ii,iii,iiii,iiiii,iiiiii)]-
		 	      Ua[s(i,ii,iii,iiii,iiiii,iiiiii)]);
                if(err>maxerr) maxerr = err;
                avgerr = avgerr + err;
                System.out.println("ug["+i+","+ii+","+iii+","+iiii+","+
                                   iiiii+","+iiiiii+"]="+
                                   U[s(i,ii,iii,iiii,iiiii,iiiiii)]+
		                   ", Ua="+Ua[s(i,ii,iii,iiii,iiiii,iiiiii)]+
                                   ", err="+err);
              } // iiiiii 
            } // iiiii 
          } // 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)nvwxyzt));
    System.out.println(" ");
    System.out.println("end pde46h_eq.java ");
  } // end pde46h_eq

  // utility functions
  int s(int i, int ii, int iii, int iiii, int iiiii, int iiiiii)
       // for v,w,x,y,z,t  in loops  for(int i=0; i<nv; i++) i!=0, i!=nv-1 
       // no boundary
  {
    int k;

    k = (i-1)*(nw-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2) +
        (ii-1)*(nx-2)*(ny-2)*(nz-2)*(nt-2) +
        (iii-1)*(ny-2)*(nz-2)*(nt-2) + (iiii-1)*(nz-2)*(nt-2) +
        (iiiii-1)*(nt-2) + (iiiiii-1);
    if(k<0 || k>=nvwxyzt) System.out.println("bad s("+i+","+ii+","+iii+
				  ","+iiii+","+iiiii+","+iiiiii+") ="+k);  
    return k;
  } // end s 

  int sk(int k, int i, int ii, int iii, int iiii, int iiiii, int iiiiii)
                                                         // for v,w,x,y,z,t 
  {
    return k*(nvwxyzt+1) + (i-1)*(nw-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2) +
           (ii-1)*(nx-2)*(ny-2)*(nz-2)*(nt-2) +
           (iii-1)*(ny-2)*(nz-2)*(nt-2) + (iiii-1)*(nz-2)*(nt-2) +
           (iiiii-1)*(nt-2) + (iiiiii-1);
  } // end sk 

  int ss(int i, int ii, int iii, int iiii, int iiiii, int iiiiii,
         int j, int jj, int jjj, int jjjj, int jjjjj, int jjjjjj)
  {
    return s(i,ii,iii,iiii,iiiii,iiiiii)*(nvwxyzt+1) + 
      s(j,jj,jjj,jjjj,jjjjj,jjjjjj);
  } // end ss 


 double uana(double v, double w, double x, double y, double z, double t)
                                        //analytic solution and boundaries
  {
    return Math.sin(v+w+x+y+z+t);
  } // end uana

  double eval_derivative(int vord, int word, int xord, int yord, int zord,
                         int tord, int i, int ii, int iii, int iiii,
                         int iiiii, int iiiiii, double u[])
  {
    int j, jj, jjj, jjjj, jjjjj, jjjjjj;
    double cv[] = new double[nv];
    double cw[] = new double[nw];
    double cx[] = new double[nx];
    double cy[] = new double[ny];
    double cz[] = new double[nz];
    double ct[] = new double[nt];
    double discrete;
    double v, w, x, y, z, t;

    new rderiv(vord, nv, i, hv, cv);
    v = vg[i];
    new rderiv(word, nw, ii, hw, cw);
    w = wg[ii];
    new rderiv(xord, nx, iii, hx, cx);
    x = xg[iii];
    new rderiv(yord, ny, iiii, hy, cy);
    y = yg[iiii];
    new rderiv(zord, nz, iiiii, hz, cz);
    z = zg[iiiii];
    new rderiv(tord, nt, iiiiii, ht, ct);
    t = tg[iiiiii];
    discrete = 0.0;


    // six cases of one partial v, w, x, y, z, t to any order 
    if(vord!=0 && word==0 && xord==0 && yord==0 && zord==0 && tord==0)
      for(j=0; j<nv; j++)
        if(j==0 || j==nv-1) discrete += cv[j]*uana(vg[j],w,x,y,z,t);
        else discrete += cv[j]*u[s(j,ii,iii,iiii,iiiii,iiiiii)];
    else if(vord==0 && word!=0 && xord==0 && yord==0 && zord==0 && tord==0)
      for(jj=0; jj<nw; jj++)
        if(jj==0 || jj==nw-1) discrete += cw[jj]*uana(v,wg[jj],x,y,z,t);
        else discrete += cw[jj]*u[s(i,jj,iii,iiii,iiiii,iiiiii)];
    else if(vord==0 && word==0 && xord!=0 && yord==0 && zord==0 && tord==0)
      for(jjj=0; jjj<nx; jjj++)
        if(jjj==0 || jjj==nx-1) discrete += cx[jjj]*uana(v,w,xg[jjj],y,z,t);
        else discrete += cx[jjj]*u[s(i,ii,jjj,iiii,iiiii,iiiiii)];
    else if(vord==0 && word==0 && xord==0 && yord!=0 && zord==0 && tord==0)
      for(jjjj=0; jjjj<ny; jjjj++)
        if(jjjj==0 || jjjj==ny-1) discrete += cy[jjjj]*uana(v,w,x,yg[jjjj],z,t);
        else discrete += cy[jjjj]*u[s(i,ii,iii,jjjj,iiiii,iiiiii)];
    else if(vord==0 && word==0 && xord==0 && yord==0 && zord!=0 && tord==0)
      for(jjjjj=0; jjjjj<nz; jjjjj++)
        if(jjjjj==0 || jjjjj==nz-1) discrete += cz[jjjjj]*
                                                uana(v,w,x,y,zg[jjjjj],t);
        else discrete += cz[jjjjj]*u[s(i,ii,iii,iiii,jjjjj,iiiiii)];
    else if(vord==0 && word==0 && xord==0 && yord==0 && zord==0 && tord!=0)
      for(jjjjjj=0; jjjjjj<nt; jjjjjj++)
        if(jjjjjj==0 || jjjjjj==nt-1) discrete += ct[jjjjjj]*
                                                  uana(v,w,x,y,z,tg[jjjjjj]);
        else discrete += ct[jjjjjj]*u[s(i,ii,iii,iiii,iiiii,jjjjjj)];
    return discrete;
  } // end eval_deriv 

  void check_soln(double u[])
  {
    // check uvvvv(v,w,x,y,z,t) + uwwww(v,w,x,y,z,t) + uxxxx(v,w,x,y,z,t) +
    //       uyyyy(v,w,x,y,z,t) + uzzzz(v,w,x,y,z,t) + utttt(v,w,x,y,z,t) +
    //       2 uvv(v,w,x,y,z,t) + 2 uww(v,w,x,y,z,t) + 2 uxx(v,w,x,y,z,t) +
    //       2 uyy(v,w,x,y,z,t) + 2 uzz(v,w,x,y,z,t) + 2 utt(v,w,x,y,z,t) +
    //       6 u(v,w,x,y,z,t) = 0
  
    double val, err, smaxerr;
    int i, ii, iii, iiii, iiiii, iiiiii;
    int ib=0, iib=0, iiib=0, iiiib=0, iiiiib=0, iiiiiib=0;

    smaxerr = 0.0;
    for(i=1; i<nv-1; i++)
    {
      for(ii=1; ii<nw-1; ii++)
      {
        for(iii=1; iii<nx-1; iii++)
        {
          for(iiii=1; iiii<ny-1; iiii++)
          {
            for(iiiii=1; iiiii<nz-1; iiiii++)
            {
              for(iiiiii=1; iiiiii<nt-1; iiiiii++)
              {
      val = 0.0;
      // uvvvv(v,w,x,y,z,t) 
      val += eval_derivative(4,0,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      // + uwwww(v,w,x,y,z,t) 
      val += eval_derivative(0,4,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      // uxxxx(v,w,x,y,z,t) 
      val += eval_derivative(0,0,4,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      // + uyyyy(v,w,x,y,z,t) 
      val += eval_derivative(0,0,0,4,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      // + uzzzz(v,w,x,y,z,t) 
      val += eval_derivative(0,0,0,0,4,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      // + utttt(v,w,x,y,z,t) 
      val += eval_derivative(0,0,0,0,0,4,i,ii,iii,iiii,iiiii,iiiiii,u);
      // + 2 uvv(v,w,x,y,z,t) 
      val += 2.0*eval_derivative(2,0,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      // + 2 uww(v,w,x,y,z,t) 
      val += 2.0*eval_derivative(0,2,0,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      // + 2 uxx(v,w,x,y,z,t) 
      val += 2.0*eval_derivative(0,0,2,0,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      // + 2 uyy(v,w,x,y,z,t) 
      val += 2.0*eval_derivative(0,0,0,2,0,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      // + 2 uzz(v,w,x,y,z,t) 
      val += 2.0*eval_derivative(0,0,0,0,2,0,i,ii,iii,iiii,iiiii,iiiiii,u);
      // + 2 utt(v,w,x,y,z,t) 
      val += 2.0*eval_derivative(0,0,0,0,0,2,i,ii,iii,iiii,iiiii,iiiiii,u);
      // + 6 u(v,w,x,y,z,t) 
      val += 6.0*uana(vg[i],wg[ii],xg[iii],yg[iiii],zg[iiiii],tg[iiiiii]);
      // -       F(v,w,x,y,z,t)  zero 
      err = Math.abs(val);
      if(err>smaxerr) 
      {
        smaxerr = err;
        ib = i;
        iib = ii;
        iiib = iii;
        iiiib = iiii;
        iiiiib = iiiii;
        iiiiiib = iiiiii;
      }
              } // end iiiii 
            } // end iiiii 
          } // end iiii 
        } // end iii 
      } // end ii 
    } // end i 
    System.out.println("check_soln maxerr="+smaxerr);
    System.out.println("at("+ib+","+iib+","+iiib+","+iiiib+
                       ","+iiiiib+","+iiiiiib+")");
  } // end check_soln 

  void write_soln(double u[]) // for plot6d_gl 
  {
    double v, w, x, y, z, t;
    try
    {
      FileWriter fout = new FileWriter("pde46h_eq_java.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde46h_eq_java.dat");
      fileout.write("6 "+nv+" "+nw+" "+nx+" "+ny+" "+nz+" "+nt+"\n");
      for(int i=0; i<nv; i++)
      {
	v = vg[i];
        for(int ii=0; ii<nw; ii++)
        {
          w = wg[ii];
          for(int iii=0; iii<nx; iii++)
          {
            x = xg[iii];
            for(int iiii=0; iiii<ny; iiii++)
            {
              y = yg[iiii];
              for(int iiiii=0; iiiii<nz; iiiii++)
              {
                z = zg[iiiii];
                for(int iiiiii=0; iiiiii<nt; iiiiii++)
                {
                  t = tg[iiiiii];
		  if(i==0 || i==nv-1 || ii==0 || ii==nw-1 ||
                     iii==0 || iii==nx-1 || iiii==0 || iiii==ny-1 ||
                     iiiii==0 || iiiii==nz-1 || iiiiii==0 || iiiiii==nt-1)
                    fileout.write(" "+v+" "+w+" "+x+" "+y+" "+z+
		  		  " "+t+" "+uana(v,w,x,y,z,t)+" \n");
                  else
                    fileout.write(" "+v+" "+w+" "+x+" "+y+" "+z+
                                  " "+t+" "+u[s(i,ii,iii,iiii,iiiii,iiiiii)]+
                                  " \n");
	        } // end iiiiii
	      } // end iiiii
	    } // end iiii
	  } // end iii
        } // end ii
      } // end i

      fileout.close();
      System.out.println("finished writing pde46h_eq_java.dat");
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("pde46h_eq_java.dat  not opened");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("pde46h_eq_java.java output-file error");
    } // end try
  } // end write_soln
 


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