abstract class f_xy extends Mathd
{
  // single function  multi independent variable
  // a single value is returned indiced to equation_ref
abstract public double func(double x,double y);
public double func(double x[])
{return func(x[0],x[1]);}
public double[] dfunc(double x1,double y1)
{ double x[]={x1,y1};
	int n=x.length;
 double c[]=new double[n];
 for(int i=0;i<n;i++)
 {c[i]=dfunc(x,i);}
 return c;
}

public double[][] d2func(double x1,double y1)
{double x[]={x1,y1};
	int n=x.length;
 double c[][]=new double[n][n];

 for(int i=0;i<n;i++)
 {for(int j=0;j<n;j++)
   {   int xref[]={i,j};
	   c[i][j]=d2func(x,xref);
   }
 }
 return c;
}

public double d2func(double x[],int x_ref[])
{// derivative of the function with respect to x_ref
double h0=0.0256;
int i,m;
int n=7;
double f1,f2;
double x1[];
x1=new double[x.length];
double x2[];
x2=new double[x.length];
for(i=0;i<x.length;i++)
{
x1[i]=x[i];
x2[i]=x[i];
}
//derivative of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
//vector<double> h(n,0);
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++)
{
h[i]=h0*Math.pow(r,i);
}

for(i=0;i<n;i++)
{
x1[x_ref[1]]+=h[i];
x2[x_ref[1]]-=h[i];
f1=dfunc(x1,x_ref[0]);
f2=dfunc(x2,x_ref[0]);
T[i][0]=( f1 - f2)/(2.0*h[i]);
x1[x_ref[1]]=x[x_ref[1]];
x2[x_ref[1]]=x[x_ref[1]];
}
for(m=1;m<n;m++)
{
  for(i=0;i<n-m;i++)
  {
  T[i][m]=(h[i]*h[i]*T[i+1][m-1] - h[i+m]*h[i+m]*T[i][m-1])/(h[i]*h[i]
  - h[i+m]*h[i+m]);
  }
}
double xx=T[0][n-1];
return xx;
}
public double dfunc(double x[],int x_ref)
{ // derivative of the function with respect to x_ref
double h0=0.0256;
int i,m;
int n=7;
double f1,f2;
double x1[];
x1=new double[x.length];
double x2[];
x2=new double[x.length];
for(i=0;i<x.length;i++)
{
x1[i]=x[i];
x2[i]=x[i];
}
//derivative of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
//vector<double> h(n,0);
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++)
{
h[i]=h0*Math.pow(r,i);
}

for(i=0;i<n;i++)
{
x1[x_ref]+=h[i];
x2[x_ref]-=h[i];
f1=func(x1);
f2=func(x2);
T[i][0]=( f1 - f2)/(2.0*h[i]);
x1[x_ref]=x[x_ref];
x2[x_ref]=x[x_ref];
}
for(m=1;m<n;m++)
{
  for(i=0;i<n-m;i++)
  {
  T[i][m]=(h[i]*h[i]*T[i+1][m-1] - h[i+m]*h[i+m]*T[i][m-1])/(h[i]*h[i]
  - h[i+m]*h[i+m]);
  }
}
double xx=T[0][n-1];
return xx;
}
 
}
