/* pde49hn_eq.c  homogeneous Biharmonic 4th order 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[i8], q[i9]
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include "deriv.h"
#include "simeq.h"

#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))
#undef  min
#define min(a,b) ((a)<(b)?(a):(b))
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))


#define nx  5 
#define ny  5
#define nz  5
#define nt  5
#define nu  5
#define nv  5
#define nw  5
#define np  5
#define nq  5
#define ndof (nx-2)*(ny-2)*(nz-2)*(nt-2)*(nu-2)*(nv-2)*(nw-2)*(np-2)*(nq-2)

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 */
{
  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) printf("bads ndof=%d, s(%d %d %d %d %d %d %d %d %d\n",
			    ndof, i,ii,i3,i4,i5,i6,i7,i8,i9);  
  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 */

static double A[ndof*(ndof+1)]; /* last column is RHS */
static double U[ndof];
static double xg[nx];
static double yg[ny];
static double zg[nz];
static double tg[nt];
static double ug[nu];
static double vg[nv];
static double wg[nw];
static double pg[np];
static double qg[nq];
static double Ua[ndof];
static double xmin = 0.0; /* problem parameters */
static double xmax = 0.5;
static double ymin = 0.0;
static double ymax = 0.5;
static double zmin = 0.0;
static double zmax = 0.5;
static double tmin = 0.0;
static double tmax = 0.5;
static double umin = 0.0;
static double umax = 0.5;
static double vmin = 0.0;
static double vmax = 0.5;
static double wmin = 0.0;
static double wmax = 0.5;
static double pmin = 0.0;
static double pmax = 0.5;
static double qmin = 0.0;
static double qmax = 0.5;
static double hx, hy, hz, ht, hu, hv, hw, hp, hq;

static 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 sin(x+y+z+t+u+v+w+p+q);
}

static double f(double x, double y, double z, double t, double u,
                double v, double w, double p, double q)
                              /* RHS */
{
  return sin(x+y+z+t+u+v+w+p+q);
}

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[nx];
  double cy[ny];
  double cz[nz];
  double ct[nt];
  double cu[nu];
  double cv[nv];
  double cw[nw];
  double cp[np];
  double cq[nq];
  double discrete;
  double x, y, z, t, u, v, w, p, q;

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

  /* nine cases of one partial x, y, z, t, u, v, w, p, q 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 += cv[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<np; 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;

  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,0,4,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 = 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 */
  printf("check_soln maxerr=%g \n",smaxerr);
  printf("at(%1d,%1d,%1d,%1d,%1d,%1d,%1d,%1d,%1d) \n",
         ib, iib, i3b, i4b, i5b, i6b, i7b, i8b, i9b);
} /* end check_soln */

void write_soln(char tag[], double u[]) /* for plot9D_gl */
{
  FILE* outp;
  char name1[32] = "pde49hn_eq_c";
  char extn[] = ".dat";
  char * name;
  int i, ii, i3, i4, i5, i6, i7, i8, i9;

  strcat(name1,tag);
  name = strcat(name1,extn);
  printf("writing solution to %s \n", name);
  outp = fopen(name,"w");
  if(outp==NULL)
  {
    printf("ERROR unable to open %s for writing \n", name);
    return;
    exit(1);
  }
  fprintf(outp, "9 %2d %2d %2d %2d %2d %2d %2d %2d %2d\n",
                   nx, ny, nz, nt, nu, nv, nw, np, nq );
  for(i=0; i<nx; i++)
  {
    for(ii=0; ii<ny; ii++)
    {
      for(i3=0; i3<nz; i3++)
      {
        for(i4=0; i4<nt; i4++)
        {
          for(i5=0; i5<nu; i5++)
          {
            for(i6=0; i6<nv; i6++)
            {
              for(i7=0; i7<nw; i7++)
              {
                for(i8=0; i8<np; i8++)
                {
                  for(i9=0; i9<nq; 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 )
                    fprintf(outp, "%f %f %f %f %f %f %f %f %f %f\n",
		            xg[i],yg[ii],zg[i3],tg[i4],ug[i5],vg[i6],
	                    wg[i7],pg[i8],qg[i9],
                            Ub(xg[i],yg[ii],zg[i3],tg[i4],ug[i5],
			       vg[i6],wg[i7],pg[i8],qg[i9]));
                  else
                    fprintf(outp, "%f %f %f %f %f %f %f %f %f %f\n",
		            xg[i],yg[ii],zg[i3],tg[i4],ug[i5],vg[i6],
		            wg[i7],pg[i8],qg[i9],
                            U[s(i,ii,i3,i4,i5,i6,i7,i8,i9)]);
	          } // i9
	        } // i8
	      } // i7
	    } // i6
          } // i5
        } // i4
      } // i3
    } // ii
  } // i
  fclose(outp);
  printf("solution file %s written \n", name);
} /* end write_soln */
 

int main(int argc, char *argv[])
{
  double x, y, z, t, u, v, w, p, q;
  int i,ii,i3,i4,i5,i6,i7,i8,i9;
  int 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[nx], cyy[ny], czz[nz], ctt[ny], cuu[nu], cvv[nv];
  double cww[nw], cpp[np], cqq[nq]; 
  double cxxxx[nx], cyyyy[ny], czzzz[nz], ctttt[nt];
  double cuuuu[nu], cvvvv[nv], cwwww[nw], cpppp[np];
  double cqqqq[nq];
  double val, err, avgerr, maxerr;
  double time_start;
  double time_now;
  double time_last;

  printf("pde49hn_eq.c running \n");
  printf("Given Uxxxx(x,y,z,t,u,v,w,p,q) + Uyyyy(x,y,z,t,u,v,w,p,q) + \n");
  printf("      Uzzzz(x,y,z,t,u,v,w,p,q) + Utttt(x,y,z,t,u,v,w,p,q) + \n");
  printf("      Uuuuu(x,y,z,t,u,v,w,p,q) + Uvvvv(x,y,z,t,u,v,w,p,q) + \n");
  printf("      Uwwww(x,y,z,t,u,v,w,p,q) + Upppp(x,y,z,t,u,v,w,p,q) + \n");
  printf("      Uqqqq(x,y,z,t,u,v,w,p,q) + \n");
  printf("      2 Uxx(x,y,z,t,u,v,w,p,q) + 2 Uyy(x,y,z,t,u,v,w,p,q) + \n");
  printf("      2 Uzz(x,y,z,t,u,v,w,p,q) + 2 Utt(x,y,z,t,u,v,w,p,q) + \n");
  printf("      2 Uuu(x,y,z,t,u,v,w,p,q) + 2 Uvv(x,y,z,t,u,v,w,p,q) + \n");
  printf("      2 Uww(x,y,z,t,u,v,w,p,q) + 2 Upp(x,y,z,t,u,v,w,p,q) + \n");
  printf("      2 Uqq(x,y,z,t,u,v,w,p,q) + \n");
  printf("      10  u(x,y,z,t,u,v,w,p,q) = f(x,y,z,t,u,v,w,p,q) \n");
  printf("\n"); 
  printf("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax \n");
  printf("tmin<=t<=tmax umin<=u<=umax vmin<=v<=vmax \n");
  printf("wmin<=w<=wmax pmin<=p<=pmax qmin<=q<=qmax Boundaries \n");
  printf("Analytic solution U(x,y,z,t,u,v,w,p,q) =  sin(x+y+z+t+u+v+w+p+q) \n");
  printf("\n");

  printf("xmin=%g, xmax=%g, ymin=%g, ymax=%g \n", xmin, xmax, ymin, ymax);
  printf("zmin=%g, zmax=%g, tmin=%g, tmax=%g \n", zmin, zmax, tmin, tmax);
  printf("umin=%g, umax=%g, vmin=%g, vmax=%g \n", umin, umax, vmin, vmax);
  printf("wmin=%g, wmax=%g, pmin=%g, pmax=%g \n", wmin, wmax, pmin, pmax);
  printf("qmin=%g, qmax=%g \n", qmin, qmax);
  printf("nx=%d, ny=%d, nz=%d, nt=%d, nu=%d, nv=%d, nw=%d, np=%d, nq=%d \n",
         nx, ny, nz, nt, nu, nv, nw, np, nq);
  time_start = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time start = %f seconds \n", time_start);
  time_last = time_start;

  printf("x grid and analytic solution at minimums \n");
  hx = (xmax-xmin)/(double)(nx-1);
  for(i=0; i<nx; i++)
  {
    xg[i] = xmin+(double)i*hx;
    Ua[i] = Ub(xg[i],ymin,zmin,tmin,umin,vmin,wmin,pmin,qmin); /* temp */
    printf("i=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f\n",
	   i, xg[i], ymin, zmin, tmin, umin, vmin, wmin, pmin, qmin, Ua[i]);
  } /* i  x */  
  printf("\n");

  printf("y grid and analytic solution at minimums \n");
  hy = (ymax-ymin)/(double)(ny-1);
  for(ii=0; ii<ny; ii++)
  {
    yg[ii] = ymin+(double)ii*hy;
    Ua[ii] = Ub(xmin, yg[ii], zmin, tmin, umin, vmin, wmin, pmin, qmin); /* temp */
    printf("ii=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f\n",
	   ii, xmin, yg[ii], zmin, tmin, umin, vmin, wmin, pmin, qmin, Ua[ii]);
  } /* ii  y */
  printf("\n");

  printf("z grid and analytic solution at minimums \n");
  hz = (zmax-zmin)/(double)(nz-1);
  for(i3=0; i3<nz; i3++)
  {
    zg[i3] = zmin+(double)i3*hz;
    Ua[i3] = Ub(xmin, ymin, zg[i3], tmin, umin, vmin, wmin, pmin, qmin); /* temp */
    printf("i3=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f\n",
	   i3, xmin, ymin, zg[i3], tmin, umin, vmin, wmin, pmin, qmin, Ua[i3]);
  } /* i3   z */
  printf("\n");

  printf("t grid and analytic solution at minimums \n");
  ht = (tmax-tmin)/(double)(nt-1);
  for(i4=0; i4<nt; i4++)
  {
    tg[i4] = tmin+(double)i4*ht;
    Ua[i4] = Ub(xmin, ymin, zmin, tg[i4], umin, vmin, wmin, pmin, qmin); /* temp */
    printf("i4=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f\n",
	   i4, xmin, ymin, zmin, tg[i4], umin, vmin, wmin, pmin, qmin, Ua[i4]);
  } /* i4  t */
  printf("\n");

  printf("u grid and analytic solution at minimums \n");
  hu = (umax-umin)/(double)(nu-1);
  for(i5=0; i5<nu; i5++)
  {
    ug[i5] = umin+(double)i5*hu;
    Ua[i5] = Ub(xmin, ymin, zmin, tmin, ug[i5], vmin, wmin, pmin, qmin); /* temp */
    printf("i5=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f\n",
	   i5, xmin, ymin, zmin, tmin, ug[i5], vmin, wmin, pmin, qmin, Ua[i5]);
  } /* i5  u */  
  printf("\n");

  printf("v grid and analytic solution at minimums \n");
  hv = (vmax-vmin)/(nv-1);
  for(i6=0; i6<nv; i6++)
  {
    vg[i6] = vmin+(double)i6*hv;
    Ua[i6] = Ub(xmin, ymin, zmin, tmin, umin, vg[i6], wmin, pmin, qmin); /* temp */
    printf("i6=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f\n",
	   i6, xmin, ymin, zmin, tmin, umin, vg[i6], wmin, pmin, qmin, Ua[i6]);
  } /* i6  v */  
  printf("\n");

  printf("w grid and analytic solution at minimums \n");
  hw = (wmax-wmin)/(double)(nw-1);
  for(i7=0; i7<nw; i7++)
  {
    wg[i7] = wmin+(double)i7*hw;
    Ua[i7] = Ub(xmin, ymin, zmin, tmin, umin, vmin, wg[i7], pmin, qmin); /* temp */
    printf("i7=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f\n",
	   i7, xmin, ymin, zmin, tmin, umin, vmin, wg[i7], pmin, qmin, Ua[i7]);
  } /* i7  w */  
  printf("\n");

  printf("p grid and analytic solution at minimums \n");
  hp = (pmax-pmin)/(double)(np-1);
  for(i8=0; i8<np; i8++)
  {
    pg[i8] = pmin+(double)i8*hp;
    Ua[i8] = Ub(xmin, ymin, zmin, tmin, umin, vmin, wmin, pg[i8], qmin); /* temp */
    printf("i8=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f\n",
	   i8, xmin, ymin, zmin, tmin, umin, vmin, wmin, pg[i8], qmin, Ua[i8]);
  } /* i8  p */  
  printf("\n");

  printf("q grid and analytic solution at minimums \n");
  hq = (qmax-qmin)/(double)(nq-1);
  for(i9=0; i9<nq; i9++)
  {
    qg[i9] = qmin+(double)i9*hq;
    Ua[i9] = Ub(xmin, ymin, zmin, tmin, umin, vmin, wmin, pmin, qg[i9]); /* temp */
    printf("i9=%d, Ua(%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f,%6.3f)=%7.4f\n",
	   i9, xmin, ymin, zmin, tmin, umin, vmin, wmin, pmin, qg[i9], Ua[i9]);
  } /* i9  q */  
  printf("\n");


  printf("check rderiv(4,nv,0,hv,cvvvv \n");
  rderiv(4,nv,0,hv,cvvvv);
  for(i=0; i<nv; i++) printf("cvvvv[%d]=%f \n", i, cvvvv[i]);
  printf("check rderiv(4,nv,nv-1,hv,cvvvv \n");
  rderiv(4,nv,nv-1,hv,cvvvv);
  for(i=0; i<nv; i++) printf("cvvvv[%d]=%f \n", i, cvvvv[i]);
  printf("\n");

  printf("put solution in Ua vector \n");
  /* 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)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");

  printf("initialize A \n");
  /* 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)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");

  /* compute entries in A matrix, inside boundary */
  printf("compute A matrix \n");
  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 */
                  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 */
                  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 */
                  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],
		           vg[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 */
                  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 */
                  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 */
                  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 */
                  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 */
                  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 */
                  rderiv(4,nq,i9,hq,cqqqq);
                  for(j9=0; j9<nq; j9++)
                  {
                    val = cqqqq[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 */
                  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[i],yg[jj],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 */
                  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 */
                  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 */
                  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 */
                  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 */
                  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 */
                  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 */
                  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 */
                  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 */

  printf("A computed stiffness matrix \n");
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");

  printf("solve %d equations in %d unknowns \n", ndof, ndof);

  simeqb(ndof, A, U); /* A destroyed */

  printf("equations solved \n");
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("time now = %f seconds, time for previous section = %f seconds \n",
         time_now, time_now-time_last);
  time_last = time_now;
  printf("\n");
  
  printf("check PDE against known solution\n");
  check_soln(Ua);
  write_soln("1",Ua);

  printf("check PDE against computed solution\n");
  check_soln(U);
  write_soln("2",U);
  printf("\n");

  avgerr = 0.0;
  maxerr = 0.0;
  printf("          U computed, Ua analytic, error \n");
  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 = 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;
                    printf("ug[%d,%d,%d,%d,%d,%d,%d,%d,%d]=%9.4f, Ua=%9.4f, err=%e \n",
			   i,ii,i3,i4,i5,i6,i7,i8,i9,
			   U[s(i,ii,i3,i4,i5,i6,i7,i8,i9)],
			   Ua[s(i,ii,i3,i4,i5,i6,i7,i8,i9)], err);
                  } /* i9 */
                } /* i8 */
              } /* i7 */
            } /* i6 */
          } /* i5 */
        } /* i4 */
      } /* i3 */
    } /* ii */
  } /* i */
  time_now = (double)clock()/(double)CLOCKS_PER_SEC;
  printf("total CPU time  = %f seconds \n", time_now-time_start);
  printf("\n");


  printf("                                  maxerr=%e, avgerr=%e \n",
         maxerr, avgerr/(double)(ndof));
  printf("\n");
  printf("end pde49hn_eq.c \n");
  return 0;
} // end main 

// end pde49hn_eq.c