// pde23_eqn2.java accuracy check pde23_eqn.java
// 
// solve equation in pdedef232.java
// class pdedef232(), ccxx(x,y,z) ... ub(x,x,z) f(x,y,z), coefs(coef)
//
// set up and solve system of (nx-2)*(ny-2)*(nz-2) linear equations
//
// uses nuderiv.java
// uses invert.java
// uses simeq.java
// uses pdedef232.java  to define problem 

import java.text.*;
import java.io.*;

class pde23_eqn2
{
  pdedef232 D = new pdedef232();
  int nx = D.nx;
  int ny = D.ny;
  int nz = D.nz;
  int neqn = (nx-2)*(ny-2)*(nz-2); // dim-2  internal, non boundary cells
  int cs = neqn;                   // index of RHS
  double xg[] = new double[nx];    // data points
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double us[] = new double[neqn];   // solution being computed 
  double U[]  = new double[neqn];   // known solution  
  double ut[][] = new double[neqn][neqn+1]; // equations incl RHS 

pde23_eqn2()
{
  double error;         // sum of absolute exact-computed         
  double avgerror;      // average error                          
  double maxerror;      // maximum cell error                     

  double coef[] = new double[9];
  double xmin, xmax, hx;
  double ymin, ymax, hy;
  double zmin, zmax, hz;

  System.out.print("pde23_eqn2.java running \n");
  System.out.print("differential equation to solve, from pdedef232.java \n");
  System.out.print("problem defined by pdedefs232.java \n");
  System.out.print("provides D.ccxx(x,y,z), f(x,y,z), D.ub(x,y,x) coefs(coef)\n");
  System.out.print("uniform grid on cube xmin,ymin,zmin to xmax,ymax,zmax\n");
  System.out.print(" \n");

  for(int i=0; i<9; i++) coef[i] = 0.0;
  System.out.print("coef initialized \n");
  D.coefs(coef);
  System.out.print("have pde232_defs(defs) \n");
  xmin = coef[0];
  xmax = coef[1];
  // nx   = (int)coef[2];
  ymin = coef[3];
  ymax = coef[4];
  // ny   = (int)coef[5];
  zmin = coef[6];
  zmax = coef[7];
  // nz   = (int)coef[8];
  System.out.print("xmin="+xmin+", xmax="+xmax+", nx="+nx+" \n");
  System.out.print("ymin="+ymin+", ymax="+ymax+", ny="+ny+" \n");
  System.out.print("zmin="+zmin+", zmax="+zmax+", nz="+nz+" \n");

  System.out.println("build and solve "+neqn+" equations in "+neqn+
                     " unknowns \n");
  hx = (xmax-xmin)/(nx-1);
  System.out.print("hx="+hx+" \n"); 

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

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

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

  System.out.print("test D.ub(xg[1],yg[1],zg[1])="+
                   D.ub(xg[1],yg[1],zg[1])+" \n");
  System.out.print("test D.f(xg[1],yg[1],zg[1])="+
                   D.f(xg[1],yg[1],zg[1])+" \n");

  for(int i=0; i<neqn; i++)
  {
    for(int j=0; j<=neqn; j++)
	ut[i][j] = 0.0;
  }
  System.out.print("ut initialized \n");

  init_matrix(); // all in and out global
  System.out.print("matrix initialized \n");
  printB(); // for debugging
  check_rows(); // all in and out global 

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

  check_soln();
  print();

  check();

  System.out.print("finished pde23_eqn2.java \n");
} // end pde23_eqn2

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

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

void check_rows() // only for(checking 
{	// row is in ut solution coordinates 
  int row, col;
  double sum, sum1;

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

void print() // typically not known 
{
  double error = 0.0;
  double max_error = 0.0;
  double avg_error = 0.0;

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

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

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

void init_matrix()
{
  double cxx[] = new double[nx];  // nuderiv coefficients 
  double cyy[] = new double[ny];
  double czz[] = new double[nz];
  double cx[] = new double[nx]; 
  double cy[] = new double[ny];
  double cz[] = new double[nz];
  double x, y, z;
  int row;

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

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

  System.out.print("internal cells zeroed \n");
  for(int i=1; i<nx-1; i++) // inside boundary 
  {
    x = xg[i];
    for(int j=1; j<ny-1; j++)
    {
      y = yg[j];
      for(int k=1; k<nz-1; k++)
      {
	z = zg[k];
        row = S(i,j,k,ny,nz); // equation 
        // ccxx*uxx+ccyy*uyy+cczz*uzz+ccxy*uxy+ccxz*uxz+ccyz*uyz+
        // ccx*ux+ccy*uy+ccz*uz+ccc*u=D.f(x,y,z)
        // for(each term 

        // ccxx*uxx  =  d^2 D.ub(x,y,z)/dx^2 
        new nuderiv(2,nx,i,xg,cxx);
        for(int ki=0; ki<nx; ki++)
        {
          if(ki==0 || ki==nx-1)
          {
	    ut[row][cs] = ut[row][cs] - D.ccxx(x,y,z)*cxx[ki]*D.ub(xg[ki],y,z);
          }	 
          else
          {
	    ut[row][S(ki,j,k,ny,nz)] = ut[row][S(ki,j,k,ny,nz)] + 
                                       D.ccxx(x,y,z)*cxx[ki];
          }
        } // ki

        // ccyy*uyy  = d^2 D.ub(x,y,z)/dy^2 
        new nuderiv(2,ny,j,yg,cyy);
        for(int kj=0; kj<ny; kj++)
        {
          if(kj==0 || kj==ny-1)
	  {
	      ut[row][cs] = ut[row][cs] - D.ccyy(x,y,z)*cyy[kj]*D.ub(x,yg[kj],z);
          }
          else
	  {
	    ut[row][S(i,kj,k,ny,nz)] = ut[row][S(i,kj,k,ny,nz)] + 
                	               D.ccyy(x,y,z)*cyy[kj];
	  }
        } // kj

        // cczz*uzz  = d^2 D.ub(x,y,z)/dz^2 
        new nuderiv(2,nz,k,zg,czz);
        for(int kk=0; kk<nz; kk++)
        {
          if(kk==0 || kk==nz-1)
	  {
	    ut[row][cs] = ut[row][cs] - D.cczz(x,y,z)*czz[kk]*D.ub(x,y,zg[kk]);
          }
          else
	  {
	    ut[row][S(i,j,kk,ny,nz)] = ut[row][S(i,j,kk,ny,nz)] + 
		                       D.cczz(x,y,z)*czz[kk];
	  }
        } // kj

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

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

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


        // ccx*ux  =  d D.ub(x,y,z)/dx 
        new nuderiv(1,nx,i,xg,cx);
        for(int ki=0; ki<nx; ki++)
        {
          if(ki==0 || ki==nx-1)
          {
	      ut[row][cs] = ut[row][cs] - 
                            D.ccx(x,y,z)*cx[ki]*D.ub(xg[ki],y,z);
          }	 
          else
	  {
	    ut[row][S(ki,j,k,ny,nz)] = ut[row][S(ki,j,k,ny,nz)] + 
	                	       D.ccx(x,y,z)*cx[ki];
	  }
        } // ki

        // ccy*uy  = d D.ub(x,y,z)/dy 
        new nuderiv(1,ny,j,yg,cy);
        for(int kj=0; kj<ny; kj++)
        {
          if(kj==0 || kj==ny-1)
          {
	      ut[row][cs] = ut[row][cs] - D.ccy(x,y,z)*cy[kj]*D.ub(x,yg[kj],z);
          }
          else
          {
	    ut[row][S(i,kj,k,ny,nz)] = ut[row][S(i,kj,k,ny,nz)] + 
		                       D.ccy(x,y,z)*cy[kj];
          }
        } // kj

        // ccz*uz  = d D.ub(x,y,z)/dz 
        new nuderiv(1,nz,k,zg,cz);
        for(int kk=0; kk<nz; kk++)
        {
          if(kk==0 || kk==nz-1)
          {
	      ut[row][cs] = ut[row][cs] - D.ccz(x,y,z)*cz[kk]*D.ub(x,y,zg[kk]);
          }
          else
          {
	    ut[row][S(i,j,kk,ny,nz)] = ut[row][S(i,j,kk,ny,nz)] + 
		                       D.ccz(x,y,z)*cz[kk];
          }
        } // kk


        // ccc*u 
        ut[row][S(i,j,k,ny,nz)] = ut[row][S(i,j,k,ny,nz)] + D.ccc(x,y,z);


        // RHS constant 
        ut[row][cs] = ut[row][cs] + D.f(x,y,z);
        // end terms 
      } // k
    } // j 
  } // i 
} // end init_matrix 

void check_soln()  // check when solution unknown 
{                  // based on computed solution
  double cxx[] = new double[nx];  // nuderiv coefficients 
  double cyy[] = new double[ny];
  double czz[] = new double[nz];
  double cx[] = new double[nx]; 
  double cy[] = new double[ny];
  double cz[] = new double[nz];
  int row;
  int cs = (nx-2)*(ny-2)*(nz-2);
  double smaxerr, err, uxx=0.0, uyy=0.0, uzz=0.0, uxy=0.0, uxz=0.0, uyz=0.0;
  double ux=0.0, uy=0.0, uz=0.0;
  double x, y, z; // set in loops

  // check ccxx*uxx(x,y,z)+ccyy*uyy(x,y,z)+cczz*uzz(x,y,z)+
  //       ccxy*uxy(x,y,z)+ccxy*uxy(x,y,z)+ccyz*uyz(x,y,z)+
  // ccx*ux(x,y,z)+ccy*uy(x,y,z)+ccz*uz(x,y,z)+ccc*D.ub(x,y,z)=D.f(x,y,z)
  // against computed solution us[S(i,j,k,ny,nz)] at each i,j,k 
  System.out.print("check computed solution against PDE \n");
  smaxerr = 0.0;

  for(int i=1; i<nx-1; i++)  // through dof equations 
  {
    x = xg[i];
    new nuderiv(1,nx,i,xg,cx);
    new nuderiv(2,nx,i,xg,cxx);
    for(int j=1; j<ny-1; j++)
    {
      y = yg[j];
      new nuderiv(1,ny,j,yg,cy);
      new nuderiv(2,ny,j,yg,cyy);
      for(int k=1; k<nz-1; k++)
      {
	z = zg[k];
        new nuderiv(1,nz,k,zg,cz);
        new nuderiv(2,nz,k,zg,czz);

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

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

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


        uxy = 0.0;
        for(int ki=0; ki<nx; ki++) // compute numerical derivative at (i,j,k)
        {
          for(int kj=0; kj<ny; kj++)
          {
            if(ki==0 || ki==nx-1 || kj==0 || kj==ny-1)
            {
              uxy  = uxy  + cx[ki]*cy[kj]*D.ub(xg[ki],yg[kj],zg[k]); // bound
            }
            else
            {
	      uxy = uxy + cx[ki]*cy[kj]*us[S(ki,kj,k,ny,nz)];
            }
          } // kj
        } // ki

        uxz = 0.0;
        for(int ki=0; ki<nx; ki++) // compute numerical derivative at (i,j,k)
        {
          for(int kk=0; kk<nz; kk++)
          {
            if(ki==0 || ki==nx-1 || kk==0 || kk==nz-1)
            {
              uxz  = uxz  + cx[ki]*cz[kk]*D.ub(xg[ki],yg[j],zg[kk]); // bound
            }
            else
            {
              uxz = uxz + cx[ki]*cz[kk]*us[S(ki,j,kk,ny,nz)];
            }
          } // kk
        } // ki

        uyz = 0.0;
        for(int kj=0; kj<ny; kj++) // compute numerical derivative at (i,j,k) 
        {
          for(int kk=0; kk<nz; kk++)
          {
            if(kj==0 || kj==ny-1 || kk==0 || kk==nz-1)
            {
              uyz  = uyz  + cy[kj]*cz[kk]*D.ub(xg[i],yg[kj],zg[kk]); // bound
            }
            else
            {
              uyz = uyz + cy[kj]*cz[kk]*us[S(i,kj,kk,ny,nz)];
            }
          } // kk
        } // kj

        err = D.ccxx(x,y,z)*uxx + D.ccyy(x,y,z)*uyy + D.cczz(x,y,z)*uzz + 
              D.ccxy(x,y,z)*uxy + D.ccxz(x,y,z)*uxz + D.ccyz(x,y,z)*uyz +
              D.ccx(x,y,z)*ux + D.ccy(x,y,z)*uy + D.ccz(x,y,z)*uz + 
              D.ccc(x,y,z)*us[S(i,j,k,ny,nz)]
              - D.f(x,y,z); // minus right hand side 
        if(Math.abs(err) > smaxerr) smaxerr = Math.abs(err);
      } // k z
    } // j y 
  } // i x 

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

void write_soln()
{
    double x, y, z, u;
  try
  {
    FileWriter fout = new FileWriter("pde23_eqn2_java.dat");
    BufferedWriter fileout = new BufferedWriter(fout);
    System.out.println("writing pde23_eqn2_java.dat");
    fileout.write("3 "+nx+" "+ny+" "+nz+" \n");
    for(int i=0; i<nx; i++)
    {
      x = xg[i];
      for(int ii=0; ii<ny; ii++)
      {
	y = yg[ii];
	for(int iii=0; iii<nz; iii++)
	{
	  z = zg[iii];
          u = D.ub(x,y,z); // exact solution
          if(i!=0 && ii!=0 && iii!=0 && i!=nx-1 && ii!=ny-1 && iii!=nz-1)
	  {
             u =  us[S(i,ii,iii,ny,nz)];
	  }
	  fileout.write(" "+x+" "+y+" "+z+" "+u +" \n");
	} // end iii                 
      } // end ii                  
    } // end i
    fileout.close();
    System.out.println("finished writing pde23_eqn2_java.dat");
  }
  catch(FileNotFoundException exception)
  {
    System.out.println("pde23_eqn2_java.dat  not opened");
  }
  catch(Exception exception) // catches all exceptions  
  {
    System.out.println("pde23_eqn2_java.out output-file error");
  } // end try        
} // end write_soln 

  public static void main (String[] args)
  {
    new pde23_eqn2();
  }

} // end class pde23_eqn2