//======================================================
// Numerical Analysis package in java
// example to show differential equation solution
// Runge-Kutta-Verner variable step method
// Dr. Turhan Coban
// =====================================================
// dy/dx=f(x,y)
// Multiequation method with error correction
import java.io.*;
import java.util.*;

abstract class f_xi
{
  // multifunction multi independent variable
  // a single value is returned indiced to equation_ref
  // example f[0]=x[0]+sin(x[1])
  //         f[1]=x[0]*x[0]-x[1]
  // func(x,1) returns the value of f[1]
  // func(x,0) returns the value of f[0]

abstract double func(double x[],int equation_ref);
}
/*
class fm1 extends f_xi
{
//multivariable function
 double func(double x[],int x_ref)
 {
   //x[0] is x, x[1] is y
   if(x_ref==0) return -0.5*x[1];
   if(x_ref==1) return 4-0.3*x[2]-0.1*x[1];
   else return 0;
 }
}
*/
class fm1 extends f_xi
{
//multivariable function
 double func(double x[],int x_ref)
 {
   //x[0] is x, x[1] is y
   if(x_ref==0) return -0.6*x[1]+10.0*Math.exp(-(x[0]-2)*(x[0]-2)/(2*0.075*0.075));
   else return 0;
 }
}

class NA92A
{
     public static double[] doublearray(String s)
     {        
     StringTokenizer token=new StringTokenizer(s);
     int n=token.countTokens()-1;
     int m=n+1;
     double a[]=new double[m];
     int j=0;           
     while(token.hasMoreTokens())
     {
     Double ax=new Double(token.nextToken());
     a[j++]=ax.doubleValue();
     }
     return a;
     }
    
  public static double[][] RKF(f_xi fp,double x0,double xn,double f0[],double h,double eps)
  { 	  
  // Runge-Kutta Fehlberg with 6 variable step

  int M=f0.length;
  int N=0; //number of data points
  String s="";
  Vector<String> vi=new Vector<String>();
  double hi=h;
  double hh[]=new double[M];
  double deltai=1;
  double alfa=1.0;
  double delta=1;
  double deltanew=1;
  double deltamin=9.99e60;
  double deltanewmin=9.99e60;
  double hmin=h;
  double xi[]=new double[M+1];
  double xxi[]=new double[M+1];
  double k[]=new double[8];
  double fi1,fi2,fi;
  int i,j;
  double x;
  xi[0]=x0;
  x=x0;
  s=xi[0]+" ";
  for(j=0;j<M;j++)
  { xi[j+1]=f0[j];
    s+=xi[j+1]+" ";
  }
    
  while((x+h)<=xn)
  {
  xi[0]=x;	  
  for(j=1;j<=M;j++)
    {
    xxi[0]=x;
    xxi[j]=xi[j];
    k[0]=fp.func(xxi,j-1);
    xxi[0]=x+h/5.0;
    xxi[j]=xi[j]+k[0]*h/5.0;
    k[1]=fp.func(xxi,j-1);
    xxi[0]=x+h*3.0/10.0;
    xxi[j]=xi[j]+k[0]*3.0/40.0+9.0/40.0*k[1]*h;
    k[2]=fp.func(xxi,j-1);
    xxi[0]=x+3.0/5.0*h;
    xxi[j]=xi[j]+3.0/10.0*k[0]*h-9.0/10.0*k[1]*h+6.0/5.0*k[2]*h;
    k[3]=fp.func(xxi,j-1);
    xxi[0]=x+h;
    xxi[j]=xi[j]-11.0/54.0*k[0]*h+5.0/2.0*k[1]*h+70.0/27.0*k[2]*h+35.0/27.0*k[3]*h;
    k[4]=fp.func(xxi,j-1);
    xxi[0]=x+7.0/8.0*h;
    xxi[j]=xi[j]+1631.0/55296.0*k[0]*h-175.0/512.0*k[1]*h+575.0/13824.0*k[2]*h+44275.0/110592.0*k[3]*h+253.0/4096.0*k[4]*h;
    k[5]=h*fp.func(xxi,j-1);
    fi1=xi[j]+(37.0/378.0*k[0]+250.0/621.0*k[2]+125.0/594.0*k[3]+512.0/1771.0*k[5])*h;
    fi2=xi[j]+(2825.0/27648.0*k[0]+18575.0/48384.0*k[2]+13525.0/55296.0*k[3]+277.0/14336.0*k[4]+1.0/4.0*k[5])*h;
    fi=fi2;
    xi[j]=fi2;    
    delta=Math.abs(fi1-fi2);
    deltanew=Math.abs(eps*fi2);
    if(delta<=deltanew) alfa=0.2;
    else                alfa=0.25;
    hh[j-1]=h*Math.pow((deltanew/delta),alfa);
    } //end of for(j
    for(j=1;j<=M;j++)
    {if(hh[j-1]<h) {h=hh[j-1];}}
    s=xi[0]+" ";
    for(j=0;j<M;j++)
    { s+=xi[j+1]+" ";
    }
        //System.out.println(s);
    vi.add(s);
    x+=h;   
    N++;   
  } //end of while
  x-=h;
  if(x<xn)//one last additional step 
  {
  h=xn-x;
  xi[0]=xn;	  
  for(j=1;j<=M;j++)
    {
    xxi[0]=x;
    xxi[j]=xi[j];
    k[0]=fp.func(xxi,j-1);
    xxi[0]=x+h/5.0;
    xxi[j]=xi[j]+k[0]*h/5.0;
    k[1]=fp.func(xxi,j-1);
    xxi[0]=x+h*3.0/10.0;
    xxi[j]=xi[j]+k[0]*3.0/40.0+9.0/40.0*k[1]*h;
    k[2]=fp.func(xxi,j-1);
    xxi[0]=x+3.0/5.0*h;
    xxi[j]=xi[j]+3.0/10.0*k[0]*h-9.0/10.0*k[1]*h+6.0/5.0*k[2]*h;
    k[3]=fp.func(xxi,j-1);
    xxi[0]=x+h;
    xxi[j]=xi[j]-11.0/54.0*k[0]*h+5.0/2.0*k[1]*h+70.0/27.0*k[2]*h+35.0/27.0*k[3]*h;
    k[4]=fp.func(xxi,j-1);
    xxi[0]=x+7.0/8.0*h;
    xxi[j]=xi[j]+1631.0/55296.0*k[0]*h-175.0/512.0*k[1]*h+575.0/13824.0*k[2]*h+44275.0/110592.0*k[3]*h+253.0/4096.0*k[4]*h;
    k[5]=h*fp.func(xxi,j-1);
    fi1=xi[j]+(37.0/378.0*k[0]+250.0/621.0*k[2]+125.0/594.0*k[3]+512.0/1771.0*k[5])*h;
    fi2=xi[j]+(2825.0/27648.0*k[0]+18575.0/48384.0*k[2]+13525.0/55296.0*k[3]+277.0/14336.0*k[4]+1.0/4.0*k[5])*h;
    fi=fi2;
    xi[j]=fi;
    delta=Math.abs(fi1-fi2);
    deltanew=Math.abs(eps*fi2);
    if(delta<=deltanew) alfa=0.2;
    else                alfa=0.25;
    hh[j-1]=h*Math.pow((deltanew/delta),alfa);
    }
    for(j=1;j<=M;j++)
    {if(hh[j-1]<h) {h=hh[j-1];}}
    s=xi[0]+" ";
    for(j=0;j<M;j++)
    { s+=xi[j+1]+" ";
    }
    //System.out.println(s);
    vi.add(s);
    N++; 
    x+=h;     
    }  //end of if(x<xn)

   //convert vector values back to the array
       Iterator nx=vi.iterator();
       String s1="";
       
       double aa[][]=new double[M+1][N];
       i=0;
       while(nx.hasNext()) 
       {s1=(String)nx.next();
	    xi=doublearray(s1);
        for(j=0;j<=M;j++)
        { aa[j][i]=xi[j];}
        i++;
        }
   return aa;
  }  
  
  public static void main(String args[]) throws IOException
  {
  fm1 b3=new fm1();
  double y[];
  y=new double[1];
  y[0]=0.5;
  double a[][]=NA92ARKF(b3,0.0,4.0,y,0.1,1e-5);
  Text.print(Text.T(a));
  Plot pp=new Plot(a);
  //pp.setPlotType(0,24);
  pp.plot();
  }
}
