// pde45h_eq.java  homogenous Biharmonic 4th order, 4 dimensions
// solve Uwwww(w,x,y,z,t) + Uxxxx(w,x,y,z,t) + Uyyyy(w,x,y,z,t) + 
//       Uzzzz(w,x,y,z,t) + Utttt(w,x,y,z,t) + 2 Uww(w,x,y,z,t) +
//       2 Uxx(w,x,y,z,t) + 2 Uyy(w,x,y,z,t) + 2 Uzz(w,x,y,z,t) +
//       2 Utt(w,x,y,z,t) + 5 U(w,x,y,z,t) = 0
//
// boundary conditions computed using U(w,x,y,z,t)
// analytic solution is U(w,x,y,z,t) = sin(w+x+y+z+t)
//
// 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
//
 
// uses simeq.java
// uses rderiv.java
import java.text.*;
import java.io.*;

class pde45h_eq
{
  final int nw = 6; // problem parameters
  final int nx = 6;
  final int ny = 6;
  final int nz = 6;
  final int nt = 6;
  final int nwxyzt = (nw-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2);

  double A[][] = new double[nwxyzt][nwxyzt+1]; // last column is RHS 
  double U[] = new double[nwxyzt]; // computed solution
  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[nwxyzt]; // check solution
  double wmin = 0.0; // problem parameters 
  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 hw, hx, hy, hz, ht;

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

  pde45h_eq()
  {
    double w, x, y, z, t;
    int k, cs;
    // derivative point arrays, dynamic in this version 
    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 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("pde45h_eq.java running ");
    System.out.println("Given Uwwww(w,x,y,z,t) + Uxxxx(w,x,y,z,t) + Uyyyy(w,x,y,z,t) + ");
    System.out.println("      Uzzzz(w,x,y,z,t) + Utttt(w,x,y,z,t) + 2 Uww(w,x,y,z,t) + ");
    System.out.println("      2 Uxx(w,x,y,z,t) + 2 Uyy(w,x,y,z,t) + 2 Uzz(w,x,y,z,t) + ");
    System.out.println("      2 Utt(w,x,y,z,t) + 5 U(w,x,y,z,t) = 0");
    System.out.println("wmin<=w<=wmax xmin<=x<=xmax ymin<=y<=ymax Boundaries");
    System.out.println("zmin<=z<=zmax tmin<=t<=tmax Boundaries");
    System.out.println("Analytic solution U(w,x,y,z,t) =  sin(w+x+y+z+t");
    System.out.println(" ");
    System.out.println("wmin="+wmin+", wmax="+wmax+", xmin="+xmin+", xmax="+xmax);
    System.out.println("ymin="+ymin+", ymax="+ymax+", zmin="+zmin+", zmax="+zmax);
    System.out.println("tmin="+tmin+", tmax="+tmax);
    System.out.println("nw="+nw+", nx="+nx+", ny="+ny+", nz="+nz+", nt="+nt);
    time_start = (double)System.currentTimeMillis()/1000.0;
    time_last = time_start;

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

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

    System.out.println("put solution in Ua vector ");
    // put solution in Ua vector 
    for(int i=1; i<nw-1; i++)
    {
      w = wg[i];
      for(int ii=1; ii<nx-1; ii++)
      {
        x = xg[ii];
        for(int iii=1; iii<ny-1; iii++)
        {
          y = yg[iii];
          for(int iiii=1; iiii<nz-1; iiii++)
          {
            z = zg[iiii];
            for(int iiiii=1; iiiii<nt-1; iiiii++)
            {
              t = tg[iiiii];
              Ua[s(i,ii,iii,iiii,iiiii)] = uana(w,x,y,z,t);
            } // iiiii 
          } // 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<nw-1; i++)
    {
      for(int ii=1; ii<nx-1; ii++)
      {
        for(int iii=1; iii<ny-1; iii++)
        {
          for(int iiii=1; iiii<nz-1; iiii++)
          {
            for(int iiiii=1; iiiii<nt-1; iiiii++)
            {
              for(int j=1; j<nw-1; j++)
              {
                for(int jj=1; jj<nx-1; jj++)
                {
                  for(int jjj=1; jjj<ny-1; jjj++)
                  {
                    for(int jjjj=1; jjjj<nz-1; jjjj++)
                    {
                      for(int jjjjj=1; jjjjj<nt-1; jjjjj++)
                      {
		        A[s(i,ii,iii,iiii,iiiii)][s(j,jj,jjj,jjjj,jjjjj)] = 0.0;
	              } // jjjjj 
	            } // jjjj 
	          } // jjj 
	        } // jj 
              } // j 
            } // iiiii
          } // 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 = (nw-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2); // constant RHS column 
    val = 0.0;
    for(int i=1; i<nw-1; i++)
    {
      for(int ii=1; ii<nx-1; ii++)
      {
        for(int iii=1; iii<ny-1; iii++)
        {
          for(int iiii=1; iiii<nz-1; iiii++)
          {
            for(int iiiii=1; iiiii<nt-1; iiiii++)
            {
		k = s(i,ii,iii,iiii,iiiii); // row index 

              // for each term in first equation 

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

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

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

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

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

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

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

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

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

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

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

            // end terms 

	    } // end iiiii loop 
	  } // 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 "+nwxyzt+" equations in "+nwxyzt+" unknowns ");

    new simeq(nwxyzt, 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(" ");

    write_soln();

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

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

  // problem RHS, 0

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

  void check_soln(double u[])
  {
    // solve Uwwww(w,x,y,z,t) + Uxxxx(w,x,y,z,t) + Uyyyy(w,x,y,z,t) +
    //       Uzzzz(w,x,y,z,t) + Utttt(w,x,y,z,t) + 2 Uww(w,x,y,z,t) +
    //       2 Uxx(w,x,y,z,t) + 2 Uyy(w,x,y,z,t) + 2 Uzz(w,x,y,z,t) +
    //       2 Utt(w,x,y,z,t) + 5 U(w,x,y,z,t) = 0
    double val, err, smaxerr;
    int ib = 0;
    int iib = 0;
    int iiib = 0;
    int  iiiib = 0;
    int  iiiiib = 0;

    smaxerr = 0.0;
    for(int i=1; i<nw-1; i++)
    {
      for(int ii=1; ii<nx-1; ii++)
      {
        for(int iii=1; iii<ny-1; iii++)
        {
          for(int iiii=1; iiii<nz-1; iiii++)
          {
            for(int iiiii=1; iiiii<nt-1; iiiii++)
            {
              val = 0.0;
              // Uwwww(w,x,y,z,t)
              val += eval_derivative(4,0,0,0,0,i,ii,iii,iiii,iiiii,u);
              // + Uxxxx(w,x,y,z,t)
              val += eval_derivative(0,4,0,0,0,i,ii,iii,iiii,iiiii,u);
	      // + Uyyyy(w,x,y,z,t)
              val += eval_derivative(0,0,4,0,0,i,ii,iii,iiii,iiiii,u);
	      // + Uzzzz(w,x,y,z,t)
              val += eval_derivative(0,0,0,4,0,i,ii,iii,iiii,iiiii,u);
              // + Utttt(w,x,y,z,t)
              val += eval_derivative(0,0,0,0,4,i,ii,iii,iiii,iiiii,u);
	      // + 2 Uww(w,x,y,z,t)
              val += 2.0*eval_derivative(2,0,0,0,0,i,ii,iii,iiii,iiiii,u);
	      // + 2 Uxx(w,x,y,z,t)
              val += 2.0*eval_derivative(0,2,0,0,0,i,ii,iii,iiii,iiiii,u);
	      // + 2 Uyy(w,x,y,z,t)
              val += 2.0*eval_derivative(0,0,2,0,0,i,ii,iii,iiii,iiiii,u);
              // + 2 Uzz(w,x,y,z,t)
              val += 2.0*eval_derivative(0,0,0,2,0,i,ii,iii,iiii,iiiii,u);
              // + 2 Utt(w,x,y,z,t)
              val += 2.0*eval_derivative(0,0,0,0,2,i,ii,iii,iiii,iiiii,u);
	      // + 5 U(w,x,y,z,t)
              val += 5.0*uana(wg[i],xg[ii],yg[iii],zg[iiii],tg[iiiii]);
              //        F(w,x,y,z,t) should be zero
	      err = Math.abs(val);
              if(err>smaxerr)
	      {
                smaxerr = err;
                ib = i; iib = ii; iiib = iii; iiiib = iiii; iiiiib = 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);
  } // end check_soln

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

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

  double eval_derivative(int word, int xord, int yord, int zord, int tord,
                      int i, int ii, int iii, int iiii, int iiiii, double u[])
  {
    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 w, x, y, z, t;

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

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

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

//  modify for reeading nx,xg, ny,yg, nz,zg, nt,tg, bxyzt, bijkl
//  int do_input() // of  .tri  .coord  .bound
//  {
//    String fname;
//    int index, last, len;
//    // read in triangles, set ntri 
//    fname = root+".tri";
//    if(debug>2) System.out.println("about to read triangles from "+fname);
//    ntri=0;
//    minvert=999999;
//    nvert=0;
//
//    try
//    {
//      BufferedReader in = new BufferedReader(new FileReader(fname));
//      String input_line;
//      input_line = in.readLine();
//      while(input_line != null)
//      {
//	 if(debug>1) System.out.println(input_line);
//         // have a line, extract 3 integers
//         len = input_line.length();
//         if(len<5)  // not enough data, ignore blank lines
//	 {
//	   input_line = in.readLine();
//           continue;
//         }
//         index = 0;
//         while(input_line.substring(index,index+1).equals(" ")){index++;}
//         last = input_line.indexOf(' ',index+1);
//         t1[ntri] = Integer.parseInt(input_line.substring(index,last)); 
//         index = last+1;
//         while(input_line.substring(index,index+1).equals(" ")){index++;}
//         last = input_line.indexOf(' ',index+1);
//         t2[ntri] = Integer.parseInt(input_line.substring(index,last)); 
//         index = last+1;
//         while(input_line.substring(index,index+1).equals(" ")){index++;}
//         last = input_line.indexOf(' ',index+1);
//         if(last<1) last = len;
//         t3[ntri] = Integer.parseInt(input_line.substring(index,last)); 
//
//         if(debug>1) System.out.println("tri "+ntri+" has vertices "+t1[ntri]+
//                                        "  "+t2[ntri]+"  "+t3[ntri]);
//         if(t1[ntri]>nvert) nvert=t1[ntri];
//         if(t2[ntri]>nvert) nvert=t2[ntri];
//         if(t3[ntri]>nvert) nvert=t3[ntri];
//         if(t1[ntri]<minvert) minvert=t1[ntri];
//         if(t2[ntri]<minvert) minvert=t2[ntri];
//         if(t3[ntri]<minvert) minvert=t3[ntri];
//         ntri++;
//	 input_line = in.readLine();
//      }
//    }
//    catch(IOException exception)
//    {
//      System.out.println("can not open for reading "+fname);
//      System.out.println(exception);
//      return 1;
//    }
//    return 0;
//  } // end do_input 
//