// pde22a_eqp.java  2D uniform grid boundary value PDE  using parts 
//                known solution to test method         using np=7
// solve uxx + 2*uxy + 3*uyy + 4*ux + 5*uy + 6*u = f(x,y)
//       f(x,y) = 6 x^2 + 18 xy + 12 y^2 + 47 x + 62 y + 97
//       Analytic solution u(x,y)=x^2 + 2 y^2 + 3 xy + 4 x + 5 y + 6
//       inside 0,1 by 0,1 grid
//
// using 9 parts, 3 each for nx, then inside 3 each for ny
// row check in each of 9 parts
// all changes in  init_matrix, rest should be the same
//
// set up and solve system of (nx-2)*(ny-2) linear equations
//
// uses rderiv.java
// uses simeq.java
import java.text.*;
import java.io.*;

class pde22a_eqp
{
  int np = 7;           // number of points in derivative
  int nx = 20;          // includes boundary 
  int ny = 21;          // includes boundary 
  boolean debug = false;
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double cxx[] = new double[nx];    // rderiv coefficients 
  double cyy[] = new double[ny];
  double cx[] = new double[nx]; 
  double cy[] = new double[ny];
                        // nx-2 by ny-2  internal, non boundary cells 
  int neqn = (nx-2)*(ny-2);
  int cs = neqn;
  double us[] = new double[neqn];           // solution being computed 
  double ut[][] = new double[neqn][neqn+1]; // equations incl RHS 
  double xmin = 0.0;
  double xmax = 1.0;
  double ymin = 0.1;
  double ymax = 1.1;
  double hx, hy;
  double error;         // sum of absolute exact-computed         
  double avgerror;      // average error                          
  double maxerror;      // maximum cell error                     
  double time_start;
  double time_now;
  double time_last;
  DecimalFormat fe = new DecimalFormat("0.00E00");
  DecimalFormat f6 = new DecimalFormat("0.00000");

  pde22a_eqp()
  {
    System.out.println("pde22a_eqp.java running");
    System.out.println("differential equation to solve");
    System.out.println("uxx+2*uxy+3*uyy+4*ux+5*uy+6*u=f(x,y)");
    System.out.println("f(x,y) = 6 x^2 + 18 xy + 12 y^2 + 47 x + 62 y + 97");
    System.out.println("uniform grid on rectangle 0,1 to 0.1,1.1");
    System.out.println("known Solution, for testing method");
    System.out.println("u(x,y) = x^2 + 2 y^2 + 3 xy + 4 x + 5 y + 6");
    System.out.println(" ");

    time_start = (double)System.currentTimeMillis()/1000.0;
    time_last = time_start;
    hx = (xmax-xmin)/(double)(nx-1);
    hy = (ymax-ymin)/(double)(ny-1);

    for(int i=0; i<nx; i++)
    {
      xg[i] = xmin+(double)i*hx;
      System.out.println("xg("+i+")="+xg[i]);
    }
    for(int j=0; j<ny; j++)
    {
      yg[j] = ymin+(double)j*hy;
      System.out.println("yg("+j+")="+yg[j]);
    }

    System.out.println("xmin="+f6.format(xg[0])+
		       ", xmax="+f6.format(xg[nx-1])+
                       ", nx="+nx+", hx="+f6.format(hx)); 
    System.out.println("ymin="+f6.format(yg[0])+
		       ", ymax="+f6.format(yg[ny-1])+
                       ", ny="+ny+", hy="+f6.format(hy)); 
    System.out.println("u(0.5,0.5)="+
                       f6.format(u(0.5,0.5)));
    System.out.println("f(0,5,0.5)="+
                       f6.format(f(0.5,0.5)));

    init_matrix(); // big change from pde22a_eq.java
    System.out.println("matrix initialized");
    System.out.println("initial matrix, left side, upper ");
    // printB(); // lots of printout for debugging
    check_rows(); // only prints if error on row
    new simeq(neqn, ut, us);
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("time to compute solution = "+(time_now-time_last)+
                       " seconds ");
    time_last = time_now;
    check_soln();
    print();

    check();
    avgerror = error/(double)neqn;
    System.out.println("total error="+fe.format(error)+
                       ", average error="+fe.format(avgerror)+
                       ", max error="+fe.format(maxerror));
    write_soln2(us);
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("total time = "+(time_now-time_start)+
                       " seconds ");
    System.out.println("finished pde22a_eqp.java");
  } // end pde22a_eqp


  double u(double x, double y) // for boundary and optional check 
  {
    return x*x + 2.0*y*y + 3.0*x*y + 4.0*x + 5.0*y + 6.0;
  } // end u

  double f(double x, double y) // RHS of differential equation to solve 
  {
    return 6.0*x*x + 18.0*x*y + 12.0*y*y +  + 47.0*x + 62.0*y + 97.0;
  } // end f

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

  void check() // typically not known, instrumentation 
  {
    double uexact, err;
  
    error = 0.0;
    maxerror = 0.0;
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<ny-1; j++)
      { 
        uexact = u(xg[i], yg[j]);
        err = Math.abs(us[S(i,j)]-uexact);
        error = error + err;
        if(err>maxerror) maxerror=err;
      }
    }
  } // end check 

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

    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<ny-1; j++)
      {
        row = S(i,j);
        sum1 = 0.0;
        for(int ii=1; ii<nx-1; ii++)
        {
          for(int jj=1; jj<ny-1; jj++)
	  {
            col = S(ii,jj);
            sum1 = sum1 + ut[row][col]*u(xg[ii],yg[jj]); // u is check solution 
	  } // ii 
        } // jj 
        sum = sum1 - ut[row][cs];
        if(Math.abs(sum)>0.001)
        {
          System.out.println("row i="+i+", j="+j+", is bad err="+sum);
        }
      } // j 
    } //i 
  } // end check_rows 

  void print() // typically not known 
  {
    double error = 0.0;
    double max_error = 0.0;
    double avg_error = 0.0;
  
    System.out.println("exact solution u, computed us, error");
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<ny-1; j++)
      {
	error = u(xg[i],yg[j]) - us[S(i,j)];
        avg_error = avg_error + Math.abs(error);
        if(Math.abs(error)>max_error) max_error = Math.abs(error);

	System.out.println("xg["+i+"]="+f6.format(xg[i])+
	      	           ", yg["+j+"]="+f6.format(yg[j])+
			   ", u="+f6.format(u(xg[i],yg[j]))+
                           ", us="+f6.format(us[S(i,j)])+
                           ", err="+fe.format(error));
      }
    }
    System.out.println("max_error="+max_error+", avg_error="+
                       avg_error/(double)neqn);
    System.out.println(" ");
  } // end print 

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

  void init_matrix() // new, using np=7 parts  total of 9 sections
  {
    double val, sum, err;
    int row;
    int iv = 0;
    int jv = 0;
    int isec = 0;
    int jsec = 0;
    int rdi2, rdi3, rdj2, rdj3; // rdi1=i  rdj1=j
    int cs = neqn;

    System.out.println(" ");
    System.out.println("pdep using 3 part method nx="+nx);
    System.out.println("     using 3 part  ny="+ny +", neqn="+neqn+", np="+np);

    // 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 ii=1; ii<nx-1; ii++)
        {
          for(int jj=1; jj<ny-1; jj++)
	  {
	      ut[S(i,j)][S(ii,jj)] = 0.0;
          }
        }
        ut[S(i,j)][cs] = 0.0;
      }
    }

    System.out.println("internal cells zeroed ");
    for(int i=1; i<nx-1; i++) // inside boundary 
    {
      for(int j=1; j<ny-1; j++) // inside boundary
      {
        row = S(i,j); // equation 
        // uxx+2*uxy+3*uyy+4*ux+5*uy+6*u=f(x,y)
        // for(each term in PDE

        if(debug) System.out.println("row="+row+", i="+i+", j="+j);
	if(i<(np+1)/2) // first section i nx
	{
          isec = 1;
          iv = 0;
          new rderiv(2,np,i,hx,cxx);
          // uxx  =  d^2 u(x,y)/dx^2 
          for(int ki=0; ki<np; ki++)
          {
            if(ki==0)
	    {
	      ut[row][cs] -= 1.0*cxx[ki]*u(xg[0],yg[j]);
	    }	 
            else
	    {
	      ut[row][S(ki+iv,j)] += 1.0*cxx[ki];
	    }
          } // ki

          // 4*ux  = 4 d u(x,y)/dx 
          new rderiv(1,np,i,hx,cx);
          for(int ki=0; ki<np; ki++)
          {
            if(ki==0)
	    {
	      ut[row][cs] -= 4.0*cx[ki]*u(xg[0],yg[j]);
	    }	 
            else
	    {
	      ut[row][S(ki+iv,j)] += 4.0*cx[ki];
	    }
          } // ki

	  if(j<(np+1)/2) // first section j ny  in first i nx
	  {
            jsec = 1;
            jv = 0;
            // 2*uxy  =  2 d^2 u(x,y)/dxdy 
            new rderiv(1,np,i,hx,cx);
            new rderiv(1,np,j,hy,cy);
            for(int ki=0; ki<np; ki++)
            {
	      for(int kj=0; kj<np; kj++)
	      {
                if(ki==0 || kj==0)
	        {
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[ki],yg[kj]);
                }
                else
	        {
		  ut[row][S(ki+iv,kj+jv)] += 2.0*cx[ki]*cy[kj];
	        }
              } // kj 
            } // ki 

            // 3*uyy  = 3 d^2 u(x,y)/dy^2 
            new rderiv(2,np,j,hy,cyy);
            for(int kj=0; kj<np; kj++)
            {
              if(kj==0)
	      {
	        ut[row][cs] -= 3.0*cyy[kj]*u(xg[i],yg[0]);
              }
              else
	      {
	        ut[row][S(i,kj+jv)] += 3.0*cyy[kj];
	      }
            } // kj

            // 5*uy  = 5 d u(x,y)/dy 
            new rderiv(1,np,j,hy,cy);
            for(int kj=0; kj<np; kj++)
            {
              if(kj==0)
	      {
	        ut[row][cs] -= 5.0*cy[kj]*u(xg[i],yg[0]);
              }
              else
	      {
	        ut[row][S(i,kj+jv)] += 5.0*cy[kj];
	      }
            } // kj

          } // end first j ny

          if(j>=(np+1)/2 && j<ny-(np+1)/2) // second section j ny no boundary
                                           //       in first i nx lower bound
          {
            jsec = 2;
	    jv = j-(np+1)/2+1; // use S()
            rdj2 = (np-1)/2;
            // 2*uxy  =  2 d^2 u(x,y)/dxdy 
            new rderiv(1,np,i,hx,cx);
            new rderiv(1,np,rdj2,hy,cy);
            for(int ki=0; ki<np; ki++)
            {
	      for(int kj=0; kj<np; kj++)
	      {
                if(ki==0)
	        {
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[0],yg[kj+jv]);
                }
                else
	        {
		  ut[row][S(ki,kj+jv)] += 2.0*cx[ki]*cy[kj];
	        }
              } // kj 
            } // ki 

            // 3*uyy  = 3 d^2 u(x,y)/dy^2 
            new rderiv(2,np,rdj2,hy,cyy);
            for(int kj=0; kj<np; kj++)
            {
	      ut[row][S(i,kj+jv)] += 3.0*cyy[kj];
            } // kj

            // 5*uy  = 5 d u(x,y)/dy 
            new rderiv(1,np,rdj2,hy,cy);
            for(int kj=0; kj<np; kj++)
            {
	      ut[row][S(i,kj+jv)] += 5.0*cy[kj];
            } // kj
          } // end second j ny

          if(j>=ny-(np+1)/2) // third section  j ny upper boundary
                             //       in first i nx lower boundary
          {
            jsec = 3;
            jv = j-(np-(ny-j));
            rdj3 = np-(ny-j);
            // 2*uxy  =  2 d^2 u(x,y)/dxdy 
            new rderiv(1,np,i,hx,cx);
            new rderiv(1,np,rdj3,hy,cy);
            for(int ki=0; ki<np; ki++)
            {
	      for(int kj=0; kj<np; kj++)
	      {
                if(ki==0 && kj==np-1)
	        {
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[0],yg[ny-1]);
                }
                else if(ki==0)
	        {
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[0],yg[kj+jv]);
                }
                else if(kj==np-1)
	        {
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[ki+iv],yg[ny-1]);
                }
                else
	        {
		  ut[row][S(ki,kj+jv)] += 2.0*cx[ki]*cy[kj];
	        }
              } // kj 
            } // ki 

            // 3*uyy  = 3 d^2 u(x,y)/dy^2 
            new rderiv(2,np,rdj3,hy,cyy);
            for(int kj=0; kj<np; kj++)
            {
              if(kj==np-1)
	      {
	        ut[row][cs] -= 3.0*cyy[kj]*u(xg[i],yg[ny-1]);
              }
              else
	      {
	        ut[row][S(i,kj+jv)] += 3.0*cyy[kj];
	      }
            } // kj

            // 5*uy  = 5 d u(x,y)/dy 
            new rderiv(1,np,rdj3,hy,cy);
            for(int kj=0; kj<np; kj++)
            {
              if(kj==np-1)
	      {
	        ut[row][cs] -= 5.0*cy[kj]*u(xg[i],yg[ny-1]);
              }
              else
	      {
	        ut[row][S(i,kj+jv)] += 5.0*cy[kj];
	      }
            } // kj

          } // end third j ny
        } // end first i nx


        if(i>=(np+1)/2 && i<nx-(np+1)/2) // second section i nx no boundary
        {
          isec = 2;
	  iv = i-(np+1)/2 +1; // use S()
          rdi2 = (np-1)/2;
          // uxx  =  d^2 u(x,y)/dx^2 
	  new rderiv(2,np,rdi2,hx,cxx);
          for(int ki=0; ki<np; ki++) // no boundary
	  {
	    ut[row][S(ki+iv,j)] += cxx[ki];
	  }
          // 4*ux  = 4 d u(x,y)/dx 
          new rderiv(1,np,rdi2,hx,cx);
          for(int ki=0; ki<np; ki++)
	  {
	    ut[row][S(ki+iv,j)] += 4.0*cx[ki];
	  }

	  if(j<(np+1)/2) // first section j ny boundary  in second section i nx
	  {
            jsec = 1;
            // 2*uxy  =  2 d^2 u(x,y)/dxdy 
            new rderiv(1,np,rdi2,hx,cx);
            new rderiv(1,np,j,hy,cy);
            for(int ki=0; ki<np; ki++)
            {
	      for(int kj=0; kj<np; kj++)
	      {
                if(kj==0)
		{
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[ki+iv],yg[0]);
                }
                else
	        {
		  ut[row][S(ki+iv,kj)] += 2.0*cx[ki]*cy[kj];
	        }
              } // kj 
            } // ki 

            // 3*uyy  = 3 d^2 u(x,y)/dy^2 
            new rderiv(2,np,j,hy,cyy);
            for(int kj=0; kj<np; kj++)
            {
              if(kj==0)
	      {
	        ut[row][cs] -= 3.0*cyy[kj]*u(xg[i],yg[0]);
              }
              else
	      {
	        ut[row][S(i,kj)] += 3.0*cyy[kj];
	      }
            } // kj

            // 5*uy  = 5 d u(x,y)/dy 
            new rderiv(1,np,j,hy,cy);
            for(int kj=0; kj<np; kj++)
            {
              if(kj==0)
	      {
	        ut[row][cs] -= 5.0*cy[kj]*u(xg[i],yg[0]);
              }
              else
	      {
	        ut[row][S(i,kj)] += 5.0*cy[kj];
	      }
            } // kj
          } // end first j ny

          if(j>=(np+1)/2 && j<ny-(np+1)/2) // second section j ny no boundary
                                           // in second section i nx
          {
            jsec = 2;
	    jv = j-(np+1)/2+1; // use S()
            rdj2 = (np-1)/2;
            // 2*uxy  =  2 d^2 u(x,y)/dxdy 
            new rderiv(1,np,rdi2,hx,cx);
            new rderiv(1,np,rdj2,hy,cy);
            for(int ki=0; ki<np; ki++)
            {
	      for(int kj=0; kj<np; kj++)
	      {
		ut[row][S(ki+iv,kj+jv)] += 2.0*cx[ki]*cy[kj];
              } // kj 
            } // ki 

            // 3*uyy  = 3 d^2 u(x,y)/dy^2 
            new rderiv(2,np,rdj2,hy,cyy);
            for(int kj=0; kj<np; kj++)
            {
	      ut[row][S(i,kj+jv)] += 3.0*cyy[kj];
            } // kj

            // 5*uy  = 5 d u(x,y)/dy 
            new rderiv(1,np,rdj2,hy,cy);
            for(int kj=0; kj<np; kj++)
            {
	      ut[row][S(i,kj+jv)] += 5.0*cy[kj];
            } // kj
          } // end second j ny

          if(j>=ny-(np+1)/2) // third section  j ny upper boundary
                             //      in second i nx
          {
            jsec = 3;
            jv = j-(np-(ny-j)); // no  -1  S()
            rdj3 = np-(ny-j);
            // 2*uxy  =  2 d^2 u(x,y)/dxdy 
            new rderiv(1,np,rdi2,hx,cx);
            new rderiv(1,np,rdj3,hy,cy);
            for(int ki=0; ki<np; ki++)
            {
	      for(int kj=0; kj<np; kj++)
	      {
                if(kj==np-1)
	        {
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[ki+iv],yg[ny-1]);
                }
                else
	        {
		  ut[row][S(ki+iv,kj+jv)] += 2.0*cx[ki]*cy[kj];
	        }
              } // kj 
            } // ki 

            // 3*uyy  = 3 d^2 u(x,y)/dy^2 
            new rderiv(2,np,rdj3,hy,cyy);
            for(int kj=0; kj<np; kj++)
            {
              if(kj==np-1)
	      {
	        ut[row][cs] -= 3.0*cyy[kj]*u(xg[i],yg[ny-1]);
              }
              else
	      {
	        ut[row][S(i,kj+jv)] += 3.0*cyy[kj];
	      }
            } // kj

            // 5*uy  = 5 d u(x,y)/dy 
            new rderiv(1,np,rdj3,hy,cy);
            for(int kj=0; kj<np; kj++)
            {
              if(kj==np-1)
	      {
	        ut[row][cs] -= 5.0*cy[kj]*u(xg[i],yg[ny-1]);
              }
              else
	      {
	        ut[row][S(i,kj+jv)] += 5.0*cy[kj];
	      }
            } // kj
          } // end third j ny
        } // end second i nx

        if(i>=nx-(np+1)/2) // third section  i nx upper boundary
        {
          isec = 3;
          iv = i-(np-(nx-i));
          rdi3 = np-(nx-i);
          new rderiv(2,np,rdi3,hx,cxx);
          // uxx  =  d^2 u(x,y)/dx^2 
	  for(int ki=0; ki<np; ki++)
	  {
	    if(ki==np-1)
	    {
	      ut[row][cs] -= cxx[ki]*u(xg[nx-1],yg[j]); // upper boundary
	    }
	    else
	    {
		ut[row][S(ki+iv,j)] += cxx[ki];			      ;
	    }
	  } // end ip

          // 4*ux  = 4 d u(x,y)/dx 
          new rderiv(1,np,rdi3,hx,cx);
          for(int ki=0; ki<np; ki++)
          {
            if(ki==np-1)
	    {
	      ut[row][cs] -= 4.0*cx[ki]*u(xg[nx-1],yg[j]);
	    }	 
            else
	    {
	      ut[row][S(ki+iv,j)] += 4.0*cx[ki];
	    }
          } // ki


	  if(j<(np+1)/2) // first section j ny in third section i nx
	  {
            jsec = 1;
            jv = 0;
            // 2*uxy  =  2 d^2 u(x,y)/dxdy 
            new rderiv(1,np,rdi3,hx,cx);
            new rderiv(1,np,j,hy,cy);
            for(int ki=0; ki<np; ki++)
            {
	      for(int kj=0; kj<np; kj++)
	      {
                if(ki==np-1 && kj==0)
	        {
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[nx-1],yg[0]);
                }
                else if(ki==np-1)
		{
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[nx-1],yg[kj+jv]);
                }
                else if(kj==0)
	        {
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[ki+iv],yg[0]);
                }
                else
	        {
		  ut[row][S(ki+iv,kj)] += 2.0*cx[ki]*cy[kj];
	        }
              } // kj 
            } // ki 

            // 3*uyy  = 3 d^2 u(x,y)/dy^2 
            new rderiv(2,np,j,hy,cyy);
            for(int kj=0; kj<np; kj++)
            {
              if(kj==0)
	      {
	        ut[row][cs] -= 3.0*cyy[kj]*u(xg[i],yg[0]);
              }
              else
	      {
	        ut[row][S(i,kj)] += 3.0*cyy[kj];
	      }
            } // kj

            // 5*uy  = 5 d u(x,y)/dy 
            new rderiv(1,np,j,hy,cy);
            for(int kj=0; kj<np; kj++)
            {
              if(kj==0)
	      {
	        ut[row][cs] -= 5.0*cy[kj]*u(xg[i],yg[0]);
              }
              else
	      {
	        ut[row][S(i,kj)] += 5.0*cy[kj];
	      }
            } // kj
          } // end first j ny

          if(j>=(np+1)/2 && j<ny-(np+1)/2) // second section j ny no boundary
                                           // in third section i nx
          {
            jsec = 2;
	    jv = j-(np+1)/2+1; // use S()
            rdj2 = (np-1)/2;
            // 2*uxy  =  2 d^2 u(x,y)/dxdy 
            new rderiv(1,np,rdi3,hx,cx);
            new rderiv(1,np,rdj2,hy,cy);
            for(int ki=0; ki<np; ki++)
            {
	      for(int kj=0; kj<np; kj++)
	      {
                if(ki==np-1)
		{ 
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[nx-1],yg[kj+jv]);
                }
                else
	        {
		  ut[row][S(ki+iv,kj+jv)] += 2.0*cx[ki]*cy[kj];
	        }
              } // kj 
            } // ki 

            // 3*uyy  = 3 d^2 u(x,y)/dy^2 
            new rderiv(2,np,rdj2,hy,cyy);
            for(int kj=0; kj<np; kj++)
            {
	      ut[row][S(i,kj+jv)] += 3.0*cyy[kj];
            } // kj

            // 5*uy  = 5 d u(x,y)/dy 
            new rderiv(1,np,rdj2,hy,cy);
            for(int kj=0; kj<np; kj++)
            {
	      ut[row][S(i,kj+jv)] += 5.0*cy[kj];
            } // kj

          } // end second j ny

          if(j>=ny-(np+1)/2) // third section  j ny upper boundary
                             // in third section i nx
          {
            jsec = 3;
            jv = j-(np-(ny-j));
            rdj3 = np-(ny-j);
            // 2*uxy  =  2 d^2 u(x,y)/dxdy 
            new rderiv(1,np,rdi3,hx,cx);
            new rderiv(1,np,rdj3,hy,cy);
            for(int ki=0; ki<np; ki++)
            {
	      for(int kj=0; kj<np; kj++)
	      {
                if(ki==np-1 && kj==np-1)
	        {
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[nx-1],yg[ny-1]);
                }
                else if(ki==np-1)
	        {
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[nx-1],yg[kj+jv]);
                }
                else if(kj==np-1)
	        {
	          ut[row][cs] -= 2.0*cx[ki]*cy[kj]*u(xg[ki+iv],yg[ny-1]);
                }
                else
	        {
		  ut[row][S(ki+iv,kj+jv)] += 2.0*cx[ki]*cy[kj];
	        }
              } // kj 
            } // ki 

            // 3*uyy  = 3 d^2 u(x,y)/dy^2 
            new rderiv(2,np,rdj3,hy,cyy);
            for(int kj=0; kj<np; kj++)
            {
              if(kj==np-1)
	      {
	        ut[row][cs] -= 3.0*cyy[kj]*u(xg[i],yg[ny-1]);
              }
              else
	      {
	        ut[row][S(i,kj+jv)] += 3.0*cyy[kj];
	      }
            } // kj

            // 5*uy  = 5 d u(x,y)/dy 
            new rderiv(1,np,rdj3,hy,cy);
            for(int kj=0; kj<np; kj++)
            {
              if(kj==np-1)
	      {
	        ut[row][cs] -= 5.0*cy[kj]*u(xg[i],yg[ny-1]);
              }
              else
	      {
	        ut[row][S(i,kj+jv)] += 5.0*cy[kj];
	      }
            } // kj

          } // end third j ny

        } // end third i nx


        // 6u  all
	ut[row][S(i,j)] = ut[row][S(i,j)] + 6.0;

        // RHS constant all 
        ut[row][cs] = ut[row][cs] + f(xg[i],yg[j]);
        // end terms 

        if(debug) // debug  check 1 row
        {
          sum = 0.0;
          for(int ki=1; ki<nx-1; ki++)
            for(int kj=1; kj<ny-1; kj++)
	      sum += ut[row][S(ki,kj)]*u(xg[ki],yg[kj]);
	  err = sum - ut[row][cs];
          if(debug) System.out.println("is="+isec+", js="+jsec+", sum="+sum+
                                       ", row error="+err);
          if(debug) System.out.println(" ");
	}
      } // j 
    } // i 
  } // end init_matrix 

  void check_soln()  // check when solution unknown 
  {
    double smaxerr, err, ux=0.0, uy=0.0, uxy=0.0, uxx=0.0, uyy=0.0;

    // check uxx(x,y)+2uxy(x,y)+3uyy(x,y)+4ux(x,y)+5uy(x,y)+6u(x,y)=f(x,y)
    // against computed solution us[S(i,j)] at each i,j 
    System.out.println("check computed solution against PDE");
    smaxerr = 0.0;

    for(int i=1; i<nx-1; i++)  // through dof equations 
    {
      new rderiv(1,nx,i,hx,cx);
      new rderiv(2,nx,i,hx,cxx);
      for(int j=1; j<ny-1; j++)
      {
        new rderiv(1,ny,j,hy,cy);
        new rderiv(2,ny,j,hy,cyy);

        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]*u(xg[ki],yg[j]); // boundary 
            uxx = uxx + cxx[ki]*u(xg[ki],yg[j]); // boundary 
	  }
          else
	  {
            ux  = ux  + cx [ki]*us[S(ki,j)];
            uxx = uxx + cxx[ki]*us[S(ki,j)];
          }
        } // ki

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

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

        err = uxx + 2.0*uxy + 3.0*uyy + 4.0*ux + 5.0*uy + 6.0*us[S(i,j)]
              - f(xg[i],yg[j]); // minus right hand side 
        if(Math.abs(err) > smaxerr) smaxerr = Math.abs(err);
      } // j y 
    } // i x 

    System.out.println("check computed soln against PDE max error="+smaxerr);
    System.out.println("x="+xg[nx-2]+", y="+yg[ny-2]);
    System.out.println("uxx="+uxx+", uyy="+uyy+", ux="+ux+", uy="+uy);
    System.out.println("uxy="+uxy+", us="+us[S(nx-2,ny-2)]);
    System.out.println("f(x,y)="+f(xg[nx-2],yg[ny-2]));
    System.out.println(" ");
  } // end Check_Soln 

  // write_soln2  for gnuplot 
  void  write_soln2(double Ux[])
  {
    try
    {
      FileWriter fout = new FileWriter("pde22a_eqp_java.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde22a_eqp_java.dat");

      for(int i=0; i<nx; i++)
      {
        for(int j=0; j<ny; j++)
        {
          if(i==0 || i==nx-1 || j==0 || j==ny-1)
          {
            fileout.write(xg[i]+" "+yg[j]+" "+u(xg[i],yg[j])+"\n"); // boundary 
          }
          else
          { 
            fileout.write(xg[i]+" "+yg[j]+" "+Ux[S(i,j)]+"\n"); // computed solution 
          }
        } // j
        fileout.write("\n"); // blank line 
      } // i 
      fileout.close();
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("pde22a_eqp_java.dat  not opened");
    }
    catch(Exception exception)
    {
      System.out.println("pde22a_eqp_java.dat  not opened");
    }
    System.out.println("finished writing pde22a_eqp_java.dat");
  }  // end write_soln2 

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