public abstract  class fi_xi extends Mathd
{
  // multifunction multi independent variable
  // vector of dependent variables are returned
  // example f[0]=x[0]+sin(x[1])
  //         f[1]=x[0]*x[1]-x[1]
  // func(x) returns the value of f[0] and f[1]
  // as a two dimensional vector 
  abstract public double[] func(double x[]);

public double[] xp(double x[],double j,double m,double h)
{
//h = dx 
//j=-4..4 multiplication factor
//m derivative taken value
int k=x.length;
double xx[]=new double[k];
for(int i=0;i<k;i++)
{ 	if(i==m) xx[i]=x[i]+j*h; 
    else     xx[i]=x[i];
}
return xx;	
}
public double[][] dfunc(double x[])
{
double h=1.0e-4;
return dfunc(x,1,h);
}
public double[][] dfunc(double x[],int n)
{
double h=1.0e-4;
return dfunc(x,n,h);
}

public double[][] dfunc(double x[],int n,double h)
{
int k=x.length; //number of x values (dimension)
double df[][]=new double[k][k];	
for(int i=0;i<k;i++)
{ for(int j=0;j<k;j++)
  {df[i][j]=dfunc(x,n,i,j,h);}
}
return df;
}

public double dfunc(double x[],int n,int j,int m)
{
double h=1.0e-6;
return dfunc(x,n,j,m,h);
}
  
public double dfunc(double x[],int n,int j,int m,double h)
{ // derivative of x[m]
  //m derivative taken value
  //jth equation
  //nth derivative of the equation (n=0) function value
  //n=1 first derivative n=2 second derivative..
int k=x.length; //number of x values (dimension)
double hh=1/h;
double ff[][]=new double[9][3];
double xp[][]=new double[9][3];
double df;
int j1;
double x1=0;
x1=x[m];
    for(int jj=-4;jj<=4;jj++)
    {  j1=jj+4;
	   xp[j1]=xp(x,jj,m,h);
       ff[j1]=func(xp[j1]);
    }
if(n==0) 
{df=ff[4][j];}
else if(n==1)
{ df=(3.0*ff[0][j]-32.0*ff[1][j]+168.0*ff[2][j]-672.0*ff[3][j]+672.0*ff[5][j]-168.0*ff[6][j]+32.0*ff[7][j]-3.0*ff[8][j])/840.0*hh;}
else if(n==2)
{df=(-14350.0*ff[4][j]-9.0*ff[0][j]+128*ff[1][j]-1008*ff[2][j]+8064*ff[3][j]+8064.0*ff[5][j]-1008.0*ff[6][j]+128.0*ff[7][j]-9.0*ff[8][j])/5040.0*hh*hh;}
else if(n==3)
{df=(-7.0*ff[0][j]+72.0*ff[1][j]-338.0*ff[2][j]+488.0*ff[3][j]-488.0*ff[5][j]+338.0*ff[6][j]-72.0*ff[7][j]+7.0*ff[8][j])/240.0*hh*hh*hh;}
else if(n==4)
{df=(2730.0*ff[4][j]+7.0*ff[0][j]-96.0*ff[1][j]+676.0*ff[2][j]-1952.0*ff[3][j]-1952.0*ff[5][j]+676.0*ff[6][j]-96.0*ff[7][j]+7.0*ff[8][j])/240.0*hh*hh*hh*hh;}
else if(n==5)
{df=(ff[0][j]-9.0*ff[1][j]+26.0*ff[2][j]-29.0*ff[3][j]+29.0*ff[5][j]-26.0*ff[6][j]+9.0*ff[7][j]-ff[8][j])/6.0*hh*hh*hh*hh*hh;}
else if(n==6)
{df=(-150.0*ff[4][j]-ff[0][j]+12.0*ff[1][j]-52.0*ff[2][j]+116.0*ff[3][j]+116.0*ff[5][j]-52.0*ff[6][j]+12.0*ff[7][j]-ff[8][j])/4.0*hh*hh*hh*hh*hh*hh;}
else if(n==7)
{df=(-ff[0][j]+6.0*ff[1][j]-14.0*ff[2][j]+14.0*ff[3][j]-14.0*ff[5][j]+14.0*ff[6][j]-6.0*ff[7][j]+ff[8][j])/2.0*hh*hh*hh*hh*hh*hh*hh;}
else if(n==8)
{df=(70.0*ff[4][j]+ff[0][j]-8.0*ff[1][j]+28.0*ff[2][j]-56.0*ff[3][j]-56.0*ff[5][j]+28.0*ff[6][j]-8.0*ff[7][j]+ff[8][j])*hh*hh*hh*hh*hh*hh*hh*hh;}
else df=0;
return df;
}
  
}