// pde49hn_eq.java  Biharmonic PDE in nine dimensions
// solve Uxxxx(x,y,z,t,u,v,w,p,q) + Uyyyy(x,y,z,t,u,v,w,p,q) +
//       Uzzzz(x,y,z,t,u,v,w,p,q) + Utttt(x,y,z,t,u,v,w,p,q) +
//       Uuuuu(x,y,z,t,u,v,w,p,q) + Uvvvv(x,y,z,t,u,v,w,p,q) + 
//       Uwwww(x,y,z,t,u,v,w,p,q) + Upppp(x,y,z,t,u,v,w,p,q) +
//       Uqqqq(x,y,z,t,u,v,w,p,q) + 
//       2 Uxx(x,y,z,t,u,v,w,p,q) + 2 Uyy(x,y,z,t,u,v,w,p,q) +
//       2 Uzz(x,y,z,t,u,v,w,p,q) + 2 Utt(x,y,z,t,u,v,w,p,q) +
//       2 Uuu(x,y,z,t,u,v,w,p,q) + 2 Uvv(x,y,z,t,u,v,w,p,q) + 
//       2 Uww(x,y,z,t,u,v,w,p,q) + 2 Upp(x,y,z,t,u,v,w,p,q) +
//       2 Uqq(x,y,z,t,u,v,w,p,q) + 
//       10 U(x,y,z,t,u,v,w,p,q) = f(x,y,z,t,u,v,w,p,q)
//
// boundary conditions computed using Ub(x,y,z,t,u,v,w,p,q)
// analytic solution is U(x,y,z,t,u,v,w,p,q) = sin(x+y+z+t+u+v+w+p+q)
//
// subscripts in code use x[i], y[ii], z[i3], t[i4], u[i5],
//                        v[i6], w[i7], p[18], q[i9]
//
// replace continuous derivatives with discrete derivatives
// solve for U(x,y,z,t,u,v,w,p,q) numerical approximation U
// A * U = F, solve simultaneous equations for U
//
 
import java.text.*;
import java.io.*;
// uses rderiv.java
// uses simeq.java

class pde49hn_eq
{
  int nx = 5;
  int ny = 5;
  int nz = 5;
  int nt = 5;
  int nu = 5;
  int nv = 5;
  int nw = 5;
  int np = 5;
  int nq = 5;
  int ndof = (nx-2)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2)*(np-2)*(nq-2);
  double A[] = new double[ndof*(ndof+1)]; // last column is RHS 
  double U[] = new double[ndof]; // computed solution, no boundary
  double xg[] = new double[nx]; // coordinates
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double tg[] = new double[nt];
  double ug[] = new double[nu];
  double vg[] = new double[nv];
  double wg[] = new double[nw];
  double pg[] = new double[np];
  double qg[] = new double[nq];
  double Ua[] = new double[ndof]; // analytic solution for checking, no bou
  double xmin = 0.0; // problem parameters 
  double xmax = 0.5;
  double ymin = 0.0;
  double ymax = 0.5;
  double zmin = 0.0;
  double zmax = 0.5;
  double tmin = 0.0;
  double tmax = 0.5;
  double umin = 0.0;
  double umax = 0.5;
  double vmin = 0.0;
  double vmax = 0.5;
  double wmin = 0.0;
  double wmax = 0.5;
  double pmin = 0.0;
  double pmax = 0.5;
  double qmin = 0.0;
  double qmax = 0.5;
  double hx, hy, hz, ht, hu, hv, hw, hp, hq;

  pde49hn_eq()
  {
    double x, y, z, t, u, v, w, p, q;
    int    i,ii,i3,i4,i5,i6,i7,i8,i9,
	   j,jj,j3,j4,j5,j6,j7,j8,j9;
    int k, cs;
    // coefficients of terms in first equation 
    // derivative point arrays, dynamic in this version 
    double cxx[] = new double[nx];
    double cyy[] = new double[ny];
    double czz[] = new double[nz];
    double ctt[] = new double[nt];
    double cuu[] = new double[nu];
    double cvv[] = new double[nv]; 
    double cww[] = new double[nw]; 
    double cpp[] = new double[np]; 
    double cqq[] = new double[nq]; 
    double cxxxx[] = new double[nx];
    double cyyyy[] = new double[ny];
    double czzzz[] = new double[nz];
    double ctttt[] = new double[nt];
    double cuuuu[] = new double[nu];
    double cvvvv[] = new double[nv];
    double cwwww[] = new double[nw];
    double cpppp[] = new double[np];
    double cqqqq[] = new double[nq];
    double val, err, avgerr, maxerr;
    double time_start;
    double time_now;
    double time_last;

    System.out.println("pde49hn_eq.java running ");
    System.out.println("Given Uxxxx(x,y,z,t,u,v,w,p,q) + Uyyyy(x,y,z,t,u,v,w,p,q) + Uzzzz(x,y,z,t,u,v,w,p,q) +");
    System.out.println("      Utttt(x,y,z,t,u,v,w,p,q) + Uuuuu(x,y,z,t,u,v,w,p,q) + Uvvvv(x,y,z,t,u,v,w,p,q) +");
    System.out.println("      Uwwww(x,y,z,t,u,v,w,p,q) + Upppp(x,y,z,t,u,v,w,p,q) + Uqqqq(x,y,z,t,u,v,w,p,q) +");
    System.out.println("      2 Uxx(x,y,z,t,u,v,w,p,q) + 2 Uyy(x,y,z,t,u,v,w,p,q) + 2 Uzz(x,y,z,t,u,v,w,p,q) +");
    System.out.println("      2 Utt(x,y,z,t,u,v,w,p,q) + 2 Uuu(x,y,z,t,u,v,w,p,q) + 2 Uvv(x,y,z,t,u,v,w,p,q) +");
    System.out.println("      2 Uww(x,y,z,t,u,v,w,p,q) + 2 Upp(x,y,z,t,u,v,w,p,q) + 2 Uqq(x,y,z,t,u,v,w,p,q) +");
    System.out.println("      10  U(x,y,z,t,u,v,w,p,q) = f(x,y,z,t,u,v,w,p,q) ");
    System.out.println(" "); 
    System.out.println("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax ");
    System.out.println("tmin<=t<=tmax umin<=u<=umax vmin<=v<=vmax ");
    System.out.println("wmin<=w<=wmax pmin<=p<=pmax qmin<=q<=qmax  Boundaries ");
    System.out.println("Analytic solution U(x,y,z,t,u,v,w,p,q) =  sin(x+y+z+t+u+v+w+p+q) ");
    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("umin="+umin+", umax="+umax+
                       ", vmin="+vmin+", vmax="+vmax);
    System.out.println("wmin="+wmin+", wmax="+wmax+
                       ", pmin="+pmin+", pmax="+pmax);
    System.out.println("qmin="+qmin+", qmax="+qmax);
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz+", nt="+nt+
                     ", nu="+nu+", nv="+nv+", nw="+nw+", np="+np+
                     ", nq="+nq);
    time_start = (double)System.currentTimeMillis()/1000.0;
    System.out.println("time start="+time_start+" seconds");
    System.out.println("DOF, degrees of freedom, ndof="+ndof);

    time_last = time_start;

    System.out.println(" ");
    System.out.println("x grid and analytic solution at minimums ");
    hx = (xmax-xmin)/(double)(nx-1);
    i = 0;
    xg[i] = xmin;
    System.out.println("i="+i+", Ua("+xg[i]+","+ymin+","+zmin+
                       ","+tmin+","+umin+","+vmin+","+wmin+","+pmin+","+qmin+")="+
                       Ub(xg[i],ymin,zmin,tmin,umin,vmin,wmin,pmin,qmin)); 
    for(i=1; i<nx; i++)
    {
      xg[i] = xmin+(double)(i)*hx;
      System.out.println("i="+i+", Ua("+xg[i]+","+ymin+","+zmin+
                         ","+tmin+","+umin+","+umin+","+wmin+","+pmin+","+qmin+")="+
                         Ub(xg[i],ymin,zmin,tmin,umin,vmin,wmin,pmin,qmin)); 
    } // i  x   
    System.out.println(" ");

    System.out.println("y grid and analytic solution at minimums ");
    hy = (ymax-ymin)/(double)(ny-1);
    ii = 0;
    yg[ii] = ymin;
    System.out.println("ii="+ii+", Ua("+xmin+","+yg[ii]+","+zmin+
                     ","+tmin+","+umin+","+vmin+","+wmin+","+pmin+","+qmin+")="+
		       Ub(xmin,yg[ii],zmin,tmin,umin,vmin,wmin,pmin,qmin)); 
    for(ii=1; ii<ny; ii++)
    {
      yg[ii] = ymin+(double)(ii)*hy;
      System.out.println("ii="+ii+", Ua("+xmin+","+yg[ii]+","+zmin+
                         ","+tmin+","+umin+","+vmin+","+wmin+","+pmin+","+qmin+")="+
                         Ub(xmin,yg[ii],zmin,tmin,umin,vmin,wmin,pmin,qmin)); 
    } // ii  y   
    System.out.println(" ");

    System.out.println("z grid and analytic solution at minimums ");
    hz = (zmax-zmin)/(double)(nz-1);
    i3 = 0;
    zg[i3] = zmin;
    System.out.println("i3="+i3+", Ua("+xmin+","+ymin+","+zg[i3]+
                     ","+tmin+","+umin+","+vmin+","+wmin+","+pmin+","+qmin+")="+
		       Ub(xmin,ymin,zg[i3],tmin,umin,vmin,wmin,pmin,qmin)); 
    for(i3=1; i3<nz; i3++)
    {
      zg[i3] = zmin+(double)(i3)*hz;
      System.out.println("i3="+i3+", Ua("+xmin+","+ymin+","+zg[i3]+
                         ","+tmin+","+umin+","+vmin+","+wmin+","+pmin+","+qmin+")="+
                         Ub(xmin,ymin,zg[i3],tmin,umin,vmin,wmin,pmin,qmin)); 
    } // i3  z   
    System.out.println(" ");

    System.out.println("t grid and analytic solution at minimums ");
    ht = (tmax-tmin)/(double)(nt-1);
    i4 = 0;
    tg[i4] = tmin;
    System.out.println("i4="+i4+", Ua("+xmin+","+ymin+","+zmin+
                       ","+tg[i4]+","+umin+","+vmin+","+wmin+","+pmin+","+qmin+")="+
                       Ub(xmin,ymin,zmin,tg[i4],umin,vmin,wmin,pmin,qmin)); 
    for(i4=1; i4<nt; i4++)
    {
      tg[i4] = tmin+(double)(i4)*ht;
      System.out.println("i4="+i4+", Ua("+xmin+","+ymin+","+zmin+
                         ","+tg[i4]+","+umin+","+vmin+","+wmin+","+pmin+","+qmin+")="+
                         Ub(xmin,ymin,zmin,tg[i4],umin,vmin,wmin,pmin,qmin)); 
    } // i4  t 
    System.out.println(" ");

    System.out.println("u grid and analytic solution at minimums ");
    hu = (umax-umin)/(double)(nu-1);
    i5 = 0;
    ug[i5] = umin;
    System.out.println("i5="+i5+", Ua("+xmin+","+ymin+","+zmin+
                       ","+tmin+","+ug[i5]+","+vmin+","+wmin+","+pmin+","+qmin+")="+
                       Ub(xmin,ymin,zmin,tmin,ug[i5],vmin,wmin,pmin,qmin)); 
    for(i5=1; i5<nu; i5++)
    {
      ug[i5] = umin+(double)(i5)*hu;
      System.out.println("i5="+i5+", Ua("+xmin+","+ymin+","+zmin+
                         ","+tmin+","+ug[i5]+","+vmin+","+wmin+","+pmin+","+qmin+")="+
                         Ub(xmin,ymin,zmin,tmin,ug[i5],vmin,wmin,pmin,qmin)); 
    } // i5  u 
    System.out.println(" ");

    System.out.println("v grid and analytic solution at minimums ");
    hv = (vmax-vmin)/(double)(nv-1);
    i6 = 0;
    vg[i6] = vmin;
    System.out.println("i6="+i6+", Ua("+xmin+","+ymin+","+zmin+
                       ","+tmin+","+umin+","+vg[i6]+","+wmin+","+pmin+")="+
		       Ub(xmin,ymin,zmin,tmin,umin,vg[i6],wmin,pmin,qmin)); 
    for(i6=1; i6<nv; i6++)
    {
      vg[i6] = vmin+(double)(i6)*hv;
      System.out.println("i6="+i6+", Ua("+xmin+","+ymin+","+zmin+
                         ","+tmin+","+umin+","+vg[i6]+","+wmin+","+pmin+","+qmin+")="+
                         Ub(xmin,ymin,zmin,tmin,umin,vg[i6],wmin,pmin,qmin)); 
    } // i6  v 
    System.out.println(" ");

    System.out.println("w grid and analytic solution at minimums ");
    hw = (wmax-wmin)/(double)(nw-1);
    i7 = 0;
    wg[i7] = wmin;
    System.out.println("i7="+i7+", Ua("+xmin+","+ymin+","+zmin+
                       ","+tmin+","+umin+","+vmin+","+wg[i7]+","+qmin+")="+
		       Ub(xmin,ymin,zmin,tmin,umin,vmin,wg[i7],pmin,qmin)); 
    for(i7=1; i7<nw; i7++)
    {
      wg[i7] = wmin+(double)(i7)*hw;
      System.out.println("i7="+i7+", Ua("+xmin+","+ymin+","+zmin+
                         ","+tmin+","+umin+","+vmin+","+wg[i7]+","+pmin+","+qmin+")="+
		         Ub(xmin,ymin,zmin,tmin,umin,vmin,wg[i7],pmin,qmin)); 
    } // i7  w 
    System.out.println(" ");

    System.out.println("p grid and analytic solution at minimums ");
    hp = (pmax-pmin)/(double)(np-1);
    i8 = 0;
    pg[i8] = pmin;
    System.out.println("i8="+i8+", Ua("+xmin+","+ymin+","+zmin+
                       ","+tmin+","+umin+","+vmin+","+wmin+","+pg[i8]+","+qmin+")="+
		       Ub(xmin,ymin,zmin,tmin,umin,vmin,wmin,pg[i8],qmin)); 
    for(i8=1; i8<np; i8++)
    {
      pg[i8] = pmin+(double)(i8)*hp;
      System.out.println("i8="+i8+", Ua("+xmin+","+ymin+","+zmin+
                         ","+tmin+","+umin+","+vmin+","+wmin+","+pg[i8]+","+qmin+")="+
		         Ub(xmin,ymin,zmin,tmin,umin,vmin,wmin,pg[i8],qmin)); 
    } // i8  p 
    System.out.println(" ");

    System.out.println("q grid and analytic solution at minimums ");
    hq = (qmax-qmin)/(double)(nq-1);
    i9 = 0;
    qg[i9] = qmin;
    System.out.println("i9="+i9+", Ua("+xmin+","+ymin+","+zmin+
                       ","+tmin+","+umin+","+vmin+","+wmin+","+pmin+","+qg[i9]+")="+
		       Ub(xmin,ymin,zmin,tmin,umin,vmin,wmin,pmin,qg[i9])); 
    for(i9=1; i9<nq; i9++)
    {
      qg[i9] = qmin+(double)(i9)*hq;
      System.out.println("i9="+i9+", Ua("+xmin+","+ymin+","+zmin+
                         ","+tmin+","+umin+","+vmin+","+wmin+","+pmin+","+qg[i9]+")="+
		         Ub(xmin,ymin,zmin,tmin,umin,vmin,wmin,pmin,qg[i9])); 
    } // i9  q
    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<nx-1; i++)
  {
    x = xg[i];
    for(ii=1; ii<ny-1; ii++)
    {
      y = yg[ii];
      for(i3=1; i3<nz-1; i3++)
      {
        z = zg[i3];
        for(i4=1; i4<nt-1; i4++)
        {
          t = tg[i4];
          for(i5=1; i5<nu-1; i5++)
          {
            u = ug[i5];
            for(i6=1; i6<nv-1; i6++)
            {
              v = vg[i6];
              for(i7=1; i7<nw-1; i7++)
              {
                w = wg[i7];
                for(i8=1; i8<np-1; i8++)
                {
                  p = pg[i8];
                  for(i9=1; i9<nq-1; i9++)
                  {
                    q = qg[i9];
                    Ua[s(i,ii,i3,i4,i5,i6,i7,i8,i9)] =  Ub(x,y,z,t,u,v,w,p,q);
                  } // i9 
                } // i8 
              } // i7 
            } // i6 
          } // i5 
        } // i4 
      } // i3 
    } // 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<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      for(i3=1; i3<nz-1; i3++)
      {
        for(i4=1; i4<nt-1; i4++)
        {
          for(i5=1; i5<nu-1; i5++)
          {
            for(i6=1; i6<nv-1; i6++)
            {
              for(i7=1; i7<nw-1; i7++)
              {
                for(i8=1; i8<np-1; i8++)
                {
                  for(i9=1; i9<nq-1; i9++)
                  {
                    for(j=1; j<nx-1; j++)
                    {
                      for(jj=1; jj<ny-1; jj++)
                      {
                        for(j3=1; j3<nz-1; j3++)
                        {
                          for(j4=1; j4<nt-1; j4++)
                          {
                            for(j5=1; j5<nu-1; j5++)
                            {
                              for(j6=1; j6<nv-1; j6++)
                              {
                                for(j7=1; j7<nw-1; j7++)
                                {
                                  for(j8=1; j8<np-1; j8++)
                                  {
                                    for(j9=1; j9<nq-1; j9++)
                                    {
				      A[ss(i,ii,i3,i4,i5,i6,i7,i8,i9,
				           j,jj,j3,j4,j5,j6,j7,j8,j9)] = 0.0;
	                            } // j9 
	                          } // j8 
	                        } // j7 
	                      } // j6 
	                    } // j5 
	                  } // j4 
	                } // j3 
	              } // jj 
                    } // j 
                  } // i9 
                } // i8 
              } // i7 
            } // i6 
          } // i5 
        } // i4 
      } // i3 
    } // 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 = ndof; // constant RHS column 
  val = 0.0;
  for(i=1; i<nx-1; i++)
    {
    for(ii=1; ii<ny-1; ii++)
    {
      for(i3=1; i3<nz-1; i3++)
      {
        for(i4=1; i4<nt-1; i4++)
        {
          for(i5=1; i5<nu-1; i5++)
          {
            for(i6=1; i6<nv-1; i6++)
            { 
              for(i7=1; i7<nw-1; i7++)
              { 
                for(i8=1; i8<np-1; i8++)
                { 
                  for(i9=1; i9<nq-1; i9++)
                  { 
		    k = s(i,ii,i3,i4,i5,i6,i7,i8,i9); // row index 

                    // for each term in first equation 

	            // for d^4U/dx^4 
                    new rderiv(4,nx,i,hx,cxxxx);
                    for(j=0; j<nx; j++)
                    {
                      val = cxxxx[j];
                      if(j==0 || j==nx-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[j],yg[ii],zg[i3],tg[i4],ug[i5],
                             vg[i6],wg[i7],pg[i8],qg[i9]);
		      else A[sk(k,j,ii,i3,i4,i5,i6,i7,i8,i9)] += val;
	            }

	            // for d^4U/dy^4 
                    new rderiv(4,ny,ii,hy,cyyyy);
                    for(jj=0; jj<ny; jj++)
                    {
                      val = cyyyy[jj];
                      if(jj==0 || jj==ny-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[jj],zg[i3],tg[i4],ug[i5],
                             vg[i6],wg[i7],pg[i8],qg[i9]);
                      else A[sk(k,i,jj,i3,i4,i5,i6,i7,i8,i9)] += val;
	            }

	            // for d^4U/dz^4 
                    new rderiv(4,nz,i3,hz,czzzz);
                    for(j3=0; j3<nz; j3++)
                    {
                      val = czzzz[j3];
                      if(j3==0 || j3==nz-1)
		        A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[ii],zg[j3],tg[i4],ug[i5],
                             tg[i6],wg[i7],pg[i8],qg[i9]);
                      else A[sk(k,i,ii,j3,i4,i5,i6,i7,i8,i9)] += val;
	            }

	            // for d^4U/dt^4 
                    new rderiv(4,nt,i4,ht,ctttt);
                    for(j4=0; j4<nt; j4++)
                    {
                      val = ctttt[j4];
                      if(j4==0 || j4==nt-1)
	                A[k*(ndof+1)+cs] -=  val*
		          Ub(xg[i],yg[ii],zg[i3],tg[j4],ug[i5],
                             vg[i6],wg[i7],pg[i8],qg[i9]);
                      else A[sk(k,i,ii,i3,j4,i5,i6,i7,i8,i9)] += val;
	            }

	            // for d^4U/du^4 
                    new rderiv(4,nu,i5,hu,cuuuu);
                    for(j5=0; j5<nu; j5++)
                    {
                      val = cuuuu[j5];
                      if(j5==0 || j5==nu-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[j5],
                             vg[i6],wg[i7],pg[i8],qg[i9]);
                      else A[sk(k,i,ii,i3,i4,j5,i6,i7,i8,i9)] += val;
	            }

	            // for d^4U/dv^4 
                    new rderiv(4,nv,i6,hv,cvvvv);
                    for(j6=0; j6<nv; j6++)
                    {
                      val = cvvvv[j6];
                      if(j6==0 || j6==nv-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],
                             vg[j6],wg[i7],pg[i8],qg[i9]);
                      else A[sk(k,i,ii,i3,i4,i5,j6,i7,i8,i9)] += val;
	            }

	            // for d^4U/dw^4 
                    new rderiv(4,nw,i7,hw,cwwww);
                    for(j7=0; j7<nw; j7++)
                    {
                      val = cwwww[j7];
                      if(j7==0 || j7==nw-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],
			     vg[i6],wg[j7],pg[i8],qg[i9]);
                      else A[sk(k,i,ii,i3,i4,i5,i6,j7,i8,i9)] += val;
	            }

	            // for d^4U/dp^4 
                    new rderiv(4,np,i8,hp,cpppp);
                    for(j8=0; j8<np; j8++)
                    {
                      val = cpppp[j8];
                      if(j8==0 || j8==np-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],
			     vg[i6],wg[i7],pg[j8],qg[i9]);
                      else A[sk(k,i,ii,i3,i4,i5,i6,i7,j8,i9)] += val;
	            }

	            // for d^4U/dq^4 
                    new rderiv(4,nq,i9,hq,cqqqq);
                    for(j9=0; j9<np; j9++)
                    {
                      val = cpppp[j9];
                      if(j9==0 || j9==nq-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],
			     vg[i6],wg[i7],pg[i8],qg[j9]);
                      else A[sk(k,i,ii,i3,i4,i5,i6,i7,i8,j9)] += val;
	            }


	            // for 2 d^2U/dx^2 
                    new rderiv(2,nx,i,hx,cxx);
                    for(j=0; j<nx; j++)
                    {
                      val = 2.0*cxx[j];
                      if(j==0 || j==nx-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[j],yg[ii],zg[i3],tg[i4],ug[i5],
                             vg[i6],wg[i7],pg[i8],qg[i9]);
                      else A[sk(k,j,ii,i3,i4,i5,i6,i7,i8,i9)] += val;
	            }

	            // for 2 d^2U/dy^2 
                    new rderiv(2,ny,ii,hy,cyy);
                    for(jj=0; jj<ny; jj++)
                    {
                      val = 2.0*cyy[jj];
                      if(jj==0 || jj==ny-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[jj],zg[i3],tg[i4],ug[i5],
                             vg[i6],wg[i7],pg[i8],qg[i9]);
                      else A[sk(k,i,jj,i3,i4,i5,i6,i7,i8,i9)] += val;
	            }

	            // for 2 d^2U/dz^2 
                    new rderiv(2,nz,i3,hz,czz);
                    for(j3=0; j3<nz; j3++)
                    {
                      val = 2.0*czz[j3];
                      if(j3==0 || j3==nz-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[ii],zg[j3],tg[i4],ug[i5],
			     vg[i6],wg[i7],pg[i8],qg[i9]);
                      else A[sk(k,i,ii,j3,i4,i5,i6,i7,i8,i9)] += val;
	            }

	            // for 2 d^2U/dt^2 
                    new rderiv(2,nt,i4,ht,ctt);
                    for(j4=0; j4<nt; j4++)
                    {
                      val = 2.0*ctt[j4];
                      if(j4==0 || j4==nt-1)
	                A[k*(ndof+1)+cs] -=  val*
		          Ub(xg[i],yg[ii],zg[i3],tg[j4],ug[i5],
                             vg[i6],wg[i7],pg[i8],qg[i9]);
                      else A[sk(k,i,ii,i3,j4,i5,i6,i7,i8,i9)] += val;
	            }

	            // for 2 d^2U/du^2 
                    new rderiv(2,nu,i5,hu,cuu);
                    for(j5=0; j5<nu; j5++)
                    { 
                      val = 2.0*cuu[j5];
                      if(j5==0 || j5==nu-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[j5],
                             vg[i6],wg[i7],pg[i8],qg[i9]);
                      else A[sk(k,i,ii,i3,i4,j5,i6,i7,i8,i9)] += val;
	            }

	            // for 2 d^2U/dv^2 
                    new rderiv(2,nv,i6,hv,cvv);
                    for(j6=0; j6<nv; j6++)
                    {
                      val = 2.0*cvv[j6];
                      if(j6==0 || j6==nv-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],
                             vg[j6],wg[i7],pg[i8],qg[i9]);
                      else A[sk(k,i,ii,i3,i4,i5,j6,i7,i8,i9)] += val;
	            }

	            // for 2 d^2U/dw^2 
                    new rderiv(2,nw,i7,hw,cww);
                    for(j7=0; j7<nw; j7++)
                    {
                      val = 2.0*cww[j7];
                      if(j7==0 || j7==nw-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],
                             vg[i6],wg[j7],pg[i8],qg[i9]);
                      else A[sk(k,i,ii,i3,i4,i5,i6,j7,i8,i9)] += val;
	            }

	            // for 2 d^2U/dp^2 
                    new rderiv(2,np,i8,hp,cpp);
                    for(j8=0; j8<np; j8++)
                    {
                      val = 2.0*cpp[j8];
                      if(j8==0 || j8==np-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],
                             vg[i6],wg[i7],pg[j8],qg[i9]);
                      else A[sk(k,i,ii,i3,i4,i5,i6,i7,j8,i9)] += val;
	            }

	            // for 2 d^2U/dq^2 
                    new rderiv(2,nq,i9,hq,cqq);
                    for(j9=0; j9<nq; j9++)
                    {
                      val = 2.0*cqq[j9];
                      if(j9==0 || j9==nq-1)
	                A[k*(ndof+1)+cs] -= val*
		          Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],
                             vg[i6],wg[i7],pg[i8],qg[j9]);
                      else A[sk(k,i,ii,i3,i4,i5,i6,i7,i8,j9)] += val;
	            }

	            // for 10 u(x,y,z,t,u,v,w,p,q) 
		      A[sk(k,i,ii,i3,i4,i5,i6,i7,i8,i9)] += 10.0;

                    // for f(x,y,z,t,u,v,w,p,q) RHS 
	            A[k*(ndof+1)+cs] += f(xg[i],yg[ii],zg[i3],
		    		          tg[i4],ug[i5],vg[i6],wg[i7],pg[i8],qg[i9]);

                    // end terms 

	          } // end i9 loop 
	        } // end i8 loop 
	      } // end i7 loop 
	    } // end i6 loop 
	  } // end i5 loop 
	} // end i4 loop 
      } // end i3 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 "+ndof+" equations in "+ndof+" unknowns ");

    new simeq(ndof, 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<nx-1; i++)
    {
      for(ii=1; ii<ny-1; ii++)
      {
        for(i3=1; i3<nz-1; i3++)
        {
          for(i4=1; i4<nt-1; i4++)
          {
            for(i5=1; i5<nu-1; i5++)
            {
              for(i6=1; i6<nv-1; i6++)
              {
                for(i7=1; i7<nw-1; i7++)
                {
                  for(i8=1; i8<np-1; i8++)
                  {
                    for(i9=1; i9<nq-1; i9++)
                    {
		      err = Math.abs(U[s(i,ii,i3,i4,i5,i6,i7,i8,i9)]-
				     Ua[s(i,ii,i3,i4,i5,i6,i7,i8,i9)]);
                      if(err>maxerr) maxerr = err;
                      avgerr = avgerr + err;
                      System.out.println("ug["+i+","+ii+","+i3+","+i4+","+
                                         i5+","+i6+","+i7+","+i8+","+i9+"]="+
				         U[s(i,ii,i3,i4,i5,i6,i7,i8,i9)]+
				         ", Ua="+Ua[s(i,ii,i3,i4,i5,i6,i7,i8,i9)]+
                                         ", err="+err);
                    } // i9 
                  } // i8 
                } // i7 
              } // i6 
            } // i5 
          } // i4 
        } // i3 
      } // 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)ndof));
    System.out.println(" ");
    System.out.println("end pde49hn_eq.java ");
  } // end pde49hn_eq

  // utility functions
  int s(int i, int ii, int i3, int i4, int i5,
        int i6, int i7, int i8, int i9)
       // for x,y,z,t,u,v,w,p,q  in loops  for(int i=0; i<nx; i++) i!=0, i!=ny-1 
       // etc., no boundary
  {
    int k;

    k = (i -1)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2)*(np-2)*(nq-2) +
        (ii-1)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2)*(np-2)*(nq-2) +
        (i3-1)*(nt-2)*(nu-2)*(nv-2)*(nw-2)*(np-2)*(nq-2) + 
        (i4-1)*(nu-2)*(nv-2)*(nw-2)*(np-2)*(nq-2) +
        (i5-1)*(nv-2)*(nw-2)*(np-2)*(nq-2) + 
        (i6-1)*(nw-2)*(np-2)*(nq-2) + 
        (i7-1)*(np-2)*(nq-2) +
        (i8-1)*(nq-2) +
        (i9-1);
    if(k<0 || k>=ndof) System.out.println("bad s("+i+","+ii+","+i3+
		","+i4+","+i5+","+i6+","+i7+","+i8+","+i9+") ="+k);  
    return k;
  } // end s 

  int sk(int k, int i, int ii, int i3, int i4, int i5,
         int i6, int i7, int i8, int i9)   // for x,y,z,t,u,v,w,p,q 
  {
    return k*(ndof+1) + 
	(i -1)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2)*(np-2)*(nq-2) +
	(ii-1)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2)*(np-2)*(nq-2) +
	(i3-1)*(nt-2)*(nu-2)*(nv-2)*(nw-2)*(np-2)*(nq-2) + 
	(i4-1)*(nu-2)*(nv-2)*(nw-2)*(np-2)*(nq-2) +
	(i5-1)*(nv-2)*(nw-2)*(np-2)*(nq-2) + 
	(i6-1)*(nw-2)*(np-2)*(nq-2) + 
	(i7-1)*(np-2)*(nq-2) +
	(i8-1)*(nq-2) +
        (i9-1);
  } // end sk 

  int ss(int i, int ii, int i3, int i4, int i5, int i6, int i7, int i8, int i9,
	 int j, int jj, int j3, int j4, int j5, int j6, int j7, int j8, int j9)
  {
      return s(i,ii,i3,i4,i5,i6,i7,i8,i9)*(ndof+1) + 
	  s(j,jj,j3,j4,j5,j6,j7,j8,j9);
  } // end ss 

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

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

  double eval_derivative(int xord, int yord, int zord, int tord, int uord,
                         int vord, int word, int pord, int qord,
                         int i, int ii, int i3, int i4, int i5,
                         int i6, int i7, int i8, int i9, double U[])
  {
      int j, jj, j3, j4, j5, j6, j7, j8, j9;
    double cx[] = new double[nx];
    double cy[] = new double[ny];
    double cz[] = new double[nz];
    double ct[] = new double[nt];
    double cu[] = new double[nu];
    double cv[] = new double[nv];
    double cw[] = new double[nw];
    double cp[] = new double[np];
    double cq[] = new double[nq];
    double discrete;
    double x, y, z, t, u, v, w, p, q;

    new rderiv(xord, nx, i, hx, cx);
    x = xg[i];
    new rderiv(yord, ny, ii, hy, cy);
    y = yg[ii];
    new rderiv(word, nz, i3, hz, cz);
    z = zg[i3];
    new rderiv(tord, nt, i4, ht, ct);
    t = tg[i4];
    new rderiv(uord, nu, i5, hu, cu);
    u = ug[i5];
    new rderiv(vord, nv, i6, hv, cv);
    v = vg[i6];
    new rderiv(word, nw, i7, hw, cw);
    w = wg[i7];
    new rderiv(pord, np, i8, hp, cp);
    p = pg[i8];
    new rderiv(qord, nq, i9, hq, cq);
    q = qg[i9];
    discrete = 0.0;


    // seven cases of one partial x, y, z, t, u, v, w, p to any order 
    if(xord!=0 && yord==0 && zord==0 && tord==0 && uord==0 &&
       vord==0 && word==0 && pord==0 && qord==0)
      for(j=0; j<nx; j++)
        if(j==0 || j==nx-1) discrete += cx[j]*Ub(xg[j],y,z,t,u,v,w,p,q);
	else discrete += cx[j]*U[s(j,ii,i3,i4,i5,i6,i7,i8,i9)];
    else if(xord==0 && yord!=0 && zord==0 && tord==0 && uord==0 &&
            vord==0 && word==0 && pord==0 && qord==0)
      for(jj=0; jj<ny; jj++)
	if(jj==0 || jj==ny-1) discrete += cy[jj]*Ub(x,yg[jj],z,t,u,v,w,p,q);
	else discrete += cy[jj]*U[s(i,jj,i3,i4,i5,i6,i7,i8,i9)];
    else if(xord==0 && yord==0 && zord!=0 && tord==0 && uord==0 &&
            vord==0 && word==0 && pord==0 && qord==0)
      for(j3=0; j3<nz; j3++)
	if(j3==0 || j3==nz-1) discrete += cz[j3]*Ub(x,y,zg[j3],t,u,v,w,p,q);
	else discrete += cz[j3]*U[s(i,ii,j3,i4,i5,i6,i7,i8,i9)];
    else if(xord==0 && yord==0 && zord==0 && tord!=0 && uord==0 &&
            vord==0 && word==0 && pord==0 && qord==0)
      for(j4=0; j4<nt; j4++)
	if(j4==0 || j4==nt-1) discrete += ct[j4]*Ub(x,y,z,tg[j4],u,v,w,p,q);
	else discrete += ct[j4]*U[s(i,ii,i3,j4,i5,i6,i7,i8,i9)];
    else if(xord==0 && yord==0 && zord==0 && tord==0 && uord!=0 &&
            vord==0 && word==0 && pord==0 && qord==0)
      for(j5=0; j5<nu; j5++)
        if(j5==0 || j5==nu-1) discrete += cu[j5]*Ub(x,y,z,t,ug[j5],v,w,p,q);
        else discrete += cu[j5]*U[s(i,ii,i3,i4,j5,i6,i7,i8,i9)];
    else if(xord==0 && yord==0 && zord==0 && tord==0 && uord==0 &&
            vord!=0 && word==0 && pord==0 && qord==0)
      for(j6=0; j6<nv; j6++)
        if(j6==0 || j6==nv-1) discrete += cv[j6]* Ub(x,y,z,t,u,vg[j6],w,p,q);
        else discrete += ct[j6]*U[s(i,ii,i3,i4,i5,j6,i7,i8,i9)];
    else if(xord==0 && yord==0 && zord==0 && tord==0 && uord==0 &&
            vord==0 && word!=0 && pord==0 && qord==0)
      for(j7=0; j7<nw; j7++)
	  if(j7==0 || j7==nw-1) discrete += cw[j7]*Ub(x,y,z,t,u,v,wg[j7],p,q);
	  else discrete += cw[j7]*U[s(i,ii,i3,i4,i5,i6,j7,i8,i9)];
    else if(xord==0 && yord==0 && zord==0 && tord==0 && uord==0 &&
            vord==0 && word==0 && pord!=0 && qord==0)
      for(j8=0; j8<np; j8++)
	  if(j8==0 || j8==np-1) discrete += cp[j8]*Ub(x,y,z,t,u,v,w,pg[j8],q);
        else discrete += cp[j8]*U[s(i,ii,i3,i4,i5,i6,i7,j8,i9)];
    else if(xord==0 && yord==0 && zord==0 && tord==0 && uord==0 &&
            vord==0 && word==0 && pord==0 && qord!=0)
      for(j9=0; j9<nq; j9++)
	if(j9==0 || j9==nq-1) discrete += cq[j9]*Ub(x,y,z,t,u,v,w,p,qg[j9]);
        else discrete += cp[j9]*U[s(i,ii,i3,i4,i5,i6,i7,i8,j9)];
    return discrete;
  } // end eval_deriv 

  void check_soln(double u[])
  {
    // check Uxxxx(x,y,z,t,u,v,w,p,q) + Uyyyy(x,y,z,t,u,v,w,p,q) + 
    //       Uzzzz(x,y,z,t,u,v,w,p,q) + Utttt(x,y,z,t,u,v,w,p,q) +
    //       Uuuuu(x,y,z,t,u,v,w,p,q) + Uvvvv(x,y,z,t,u,v,w,p,q) +
    //       Uwwww(x,y,z,t,u,v,w,p,q) + Upppp(x,y,z,t,u,v,w,p,q) +
    //       Uqqqq(x,y,z,t,u,v,w,p,q) +
    //       2 Uxx(x,y,z,t,u,v,w,p,q) + 2 Uyy(x,y,z,t,u,v,w,p,q) +
    //       2 Uzz(x,y,z,t,u,v,w,p,q) + 2 Utt(x,y,z,t,u,v,w,p,q) +
    //       2 Uuu(x,y,z,t,u,v,w,p,q) + 2 Uvv(x,y,z,t,u,v,w,p,q) +
    //       2 Uww(x,y,z,t,u,v,w,p,q) + 2 Upp(x,y,z,t,u,v,w,p,q) +
    //       2 Uqq(x,y,z,t,u,v,w,p,q) +
    //       10  u(x,y,z,t,u,v,w,p,q) = f(x,y,z,t,u,v,w,p,q)
  
    double val, err, smaxerr;
    int i, ii, i3, i4, i5, i6, i7, i8, i9;
    int ib=0, iib=0, i3b=0, i4b=0, i5b=0, i6b=0, i7b=0, i8b=0, i9b=0;

    smaxerr = 0.0;
    for(i=1; i<nx-1; i++)
    {
      for(ii=1; ii<ny-1; ii++)
      {
        for(i3=1; i3<nz-1; i3++)
        {
          for(i4=1; i4<nt-1; i4++)
          {
            for(i5=1; i5<nu-1; i5++)
            {
              for(i6=1; i6<nv-1; i6++)
              {
                for(i7=1; i7<nw-1; i7++)
                {
                  for(i8=1; i8<np-1; i8++)
                  {
                    for(i9=1; i9<nq-1; i9++)
                    {
      val = 0.0;
      // Uxxxx(x,y,z,t,u,v,w,p,q) 
      val += eval_derivative(4,0,0,0,0,0,0,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + Uyyyy(x,y,z,t,u,v,w,p,q) 
      val += eval_derivative(0,4,0,0,0,0,0,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + Uzzzz(x,y,z,t,u,v,w,p,q) 
      val += eval_derivative(0,0,4,0,0,0,0,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + Utttt(x,y,z,t,u,v,w,p,q) 
      val += eval_derivative(0,0,0,4,0,0,0,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + Uuuuu(x,y,z,t,u,v,w,p,q) 
      val += eval_derivative(0,0,0,0,4,0,0,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + Uvvvv(x,y,z,t,u,v,w,p,q) 
      val += eval_derivative(0,0,0,0,4,0,0,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + Uwwww(x,y,z,t,u,v,w,p,q) 
      val += eval_derivative(0,0,0,0,0,0,4,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + Upppp(x,y,z,t,u,v,w,p,q) 
      val += eval_derivative(0,0,0,0,0,0,0,4,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + Uqqqq(x,y,z,t,u,v,w,p,q) 
      val += eval_derivative(0,0,0,0,0,0,0,0,4,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + 2 Uxx(x,y,z,t,u,v,w,p,q) 
      val += 2.0*eval_derivative(2,0,0,0,0,0,0,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + 2 Uyy(x,y,z,t,u,v,w,p,q) 
      val += 2.0*eval_derivative(0,2,0,0,0,0,0,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + 2 Uzz(x,y,z,t,u,v,w,p,q) 
      val += 2.0*eval_derivative(0,0,2,0,0,0,0,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + 2 Utt(x,y,z,t,u,v,w,p,q) 
      val += 2.0*eval_derivative(0,0,0,2,0,0,0,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + 2 Uuu(x,y,z,t,u,v,w,p,q) 
      val += 2.0*eval_derivative(0,0,0,0,2,0,0,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + 2 Uvv(x,y,z,t,u,v,w,p,q) 
      val += 2.0*eval_derivative(0,0,0,0,0,2,0,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + 2 Uww(x,y,z,t,u,v,w,p,q) 
      val += 2.0*eval_derivative(0,0,0,0,0,0,2,0,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + 2 Upp(x,y,z,t,u,v,w,p,q) 
      val += 2.0*eval_derivative(0,0,0,0,0,0,0,2,0,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + 2 Uqq(x,y,z,t,u,v,w,p,q) 
      val += 2.0*eval_derivative(0,0,0,0,0,0,0,0,2,i,ii,i3,i4,i5,i6,i7,i8,i9,U);
      // + 10 u(x,y,z,t,u,v,w,p,q) 
      val += 10.0*Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],vg[i6],wg[i7],pg[i8],qg[i9]);
      // - F(x,y,z,t,u,v,w,p,q)
      val -=  f(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],vg[i6],wg[i7],pg[i8],qg[i9]);
      err = Math.abs(val);
      if(err>smaxerr) 
      {
        smaxerr = err;
        ib = i;
        iib = ii;
        i3b = i3;
        i4b = i4;
        i5b = i5;
        i6b = i6;
        i7b = i7;
        i8b = i8;
        i9b = i9;
      }
                   } // end i9 
                 } // end i8 
               } // end i7 
             } // end i6 
            } // end i5 
          } // end i4 
        } // end i3 
      } // end ii 
    } // end i 
    System.out.println("check_soln maxerr="+smaxerr);
    System.out.println("at("+ib+","+iib+","+i3b+","+i4b+
                       ","+i5b+","+i6b+","+i7b+","+i8b+","+i9b+")");
  } // end check_soln 

  void write_soln(double U[]) // for plot9d_gl 
  {
    double x, y, z, t, u, v, w, p, q;
    try
    {
      FileWriter fout = new FileWriter("pde49hn_eq_java.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde49hn_eq_java.dat");
      fileout.write("9 "+nx+" "+ny+" "+nz+" "+nt+" "+nu+" "+nv+" "+nw+" "+np+" "+nq+"\n");
      for(int i=0; i<nx; i++)
      {
	x = xg[i];
        for(int ii=0; ii<ny; ii++)
        {
          y = yg[ii];
          for(int i3=0; i3<nz; i3++)
          {
            z = zg[i3];
            for(int i4=0; i4<nt; i4++)
            {
              t = tg[i4];
              for(int i5=0; i5<nu; i5++)
              {
                u = ug[i5];
                for(int i6=0; i6<nv; i6++)
                {
                  v = vg[i6];
                  for(int i7=0; i7<nw; i7++)
                  {
                    w = wg[i7];
                    for(int i8=0; i8<np; i8++)
                    {
                      p = pg[i8];
                      for(int i9=0; i9<nq; i9++)
                      {
                        q = qg[i9];
		        if(i ==0 || i ==nx-1 || ii==0 || ii==ny-1 ||
                           i3==0 || i3==nz-1 || i4==0 || i4==nt-1 ||
                           i5==0 || i5==nu-1 || i6==0 || i6==nv-1 ||
                           i7==0 || i7==nw-1 || i8==0 || i8==np-1 ||
                           i9==0 || i9==nq-1)
                        fileout.write(" "+x+" "+y+" "+z+" "+t+" "+u+
		   		     " "+v+" "+w+" "+p+" "+q+" "+
				      Ub(x,y,z,t,u,v,w,p,q)+" \n");
                      else
                        fileout.write(" "+x+" "+y+" "+z+" "+t+" "+u+" "+v+
		                      " "+w+" "+p+" "+q+" "+
                                      U[s(i,ii,i3,i4,i5,i6,i7,i8,i9)]+" \n");
	              } // end i9
	            } // end i8
	          } // end i7
	        } // end i6
	      } // end i5
	    } // end i4
	  } // end i3
        } // end ii
      } // end i

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

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