//======================================================
// Numerical Analysis package in java
// example to show differential equation solution
// Runge-Kutta-Fehlberg variable step method
// Dr. Turhan Coban
// =====================================================
// dy/dx=f(x,y)

import java.io.*;
import java.util.*;

abstract class f_xy
{
  //single function single independent variable
  // example f=x*x
  abstract double func(double x,double y);
}

class fxy2 extends f_xy
{
  double func(double x,double y)
  {return -0.6*y+10.0*Math.exp(-(x-2)*(x-2)/(2*0.075*0.075));}
}

class NA92H
{
     
  public static double[][] RKCK45(f_xy f,double x0,double y0,double xn,double h,double eps)
  { 	  
  // Runge-Kutta-Cash-Karp 4-5 variable-step
  // eps relative error
  int s=7;
double c[]=new double[s];
double a[][]=new double[s][s];
double bh[]=new double[s];
double b[]=new double[s];
double k[]=new double[s];
c[1] =  0;
c[2] =  0.2;
c[3] =  0.3;
c[4] =  0.6;
c[5] =  1;
c[6] = 7.0/8.0;

a[1][1] = 0;

a[2][1] =  0.2;

a[3][1] =  3.0/40.0;
a[3][2] =  9.0/40.0;

a[4][1] =  0.3;
a[4][2] =  -9.0/10.0;
a[4][3] =  6.0/5.0;

a[5][1] =  -11.0/54.0;
a[5][2] =  0.25;
a[5][3] =  -70.0/27.0;
a[5][4] =  35.0/27.0;

b[1] =  37.0/378.0;
b[2] =  250.0/621.0;
b[3] =  125.0/591.0;
b[4] =  512.0/1771.0;
b[5] =  0.0;

bh[1] =  2825.0/27648.0;
bh[2] =  18575.0/48384.0;
bh[3] =  13525.0/55296.0;
bh[4] =  277.0/14336.0;
bh[5] =  1.0/4.0;
  double alfa;
  int i=0;
  double delta=h;
  double deltayeni=eps*y0;
  double k1,k2,k3,k4,k5,k6;
  Vector<Double> xi=new Vector<Double>();
  Vector<Double> yi=new Vector<Double>();
  double x,y,yp1=0,yp2=0,xeski,yeski;
  double x1,y1;
  x=x0;
  y=y0;
  xeski=x;
  yeski=y;
  Double Xi=new Double(x);
  Double Yi=new Double(y);
  xi.add(Xi);
  yi.add(Yi);
  while((x+h)<=xn)
   {
   xeski=x;	
   yeski=y;
   yp1=y;
   yp2=y;   
   for(int ii=1;ii<s;ii++)
   {  x1=x+c[ii]*h;
      y1=y;
	  for(int j=1;j<ii;j++)
      {y1+=h*a[ii][j]*k[j];}
       k[ii]=f.func(x1,y1);
       yp1+=h*b[ii]*k[ii];
       yp2+=h*bh[ii]*k[ii];
   } 
   y=yp2;
   x=x+h;
   Xi=new Double(x);
   Yi=new Double(y);
   xi.add(Xi);
   yi.add(Yi);   
   delta=Math.abs(yp2-yp1);
   deltayeni=Math.abs(eps*h);
   if(delta<=deltayeni) alfa=0.2;
   else                 alfa=0.25;
   h=h*Math.pow((deltayeni/delta),alfa);
   i++;
   }
   if(x<xn)//son bir step daha
   {
   h=xn-x;
   for(int ii=1;ii<s;ii++)
   {  x1=x+c[ii]*h;
      y1=y;
	  for(int j=1;j<ii;j++)
      {y1+=h*a[ii][j]*k[j];}
       k[ii]=f.func(x1,y1);
       yp1+=h*b[ii]*k[ii];
       yp2+=h*bh[ii]*k[ii];
   } 
   y=yp2;
   //System.out.println("x="+x+"h="+h+"y="+yp1_6);
   x=xn;
   Xi=new Double(x);
   Yi=new Double(y);
   xi.add(Xi);
   yi.add(Yi);
   i++;
   }
   //vektöre aktarılan değerleri double matrise aktar
   int n=i;
        double aa[][]=new double[2][n+1];
        Enumeration nx=xi.elements();
        Enumeration ny=yi.elements();
        for(i=0;i<=n;i++)
        {
	    //String ss1=(String)nx.nextElement();
	    //String ss2=(String)ny.nextElement();    
	    Double ax=(Double)nx.nextElement();
	    Double ay=(Double)ny.nextElement();
        aa[0][i]=ax.doubleValue();
        aa[1][i]=ay.doubleValue();  
        }
   return aa;
  }
  
  
  public static void main(String args[]) throws IOException
  {
   double x0=0;
   double y0=0.5;
   fxy2 f1=new fxy2();
   double a[][]=RKCK45(f1,x0,y0,4.0,0.1,1e-3);
    Text.print(Text.T(a));
    Plot p1=new Plot(a);
    p1.setPlotType(0,24);
    p1.addSpline(a,8);
    //p1.addData(a);
    p1.plot();
  }
}
