//************************************************************************ 
// Solve the ODE x'=-(sin(t^3)+3*t^3*cos(t^3))*x on 0<t<3 with x(0)=1 
// within Runge-Kutta-Fehlberg method.
// Initial time step is given by h=0.003.
// Exact solution is: x(t)=exp(-t*sin(t^3)).
//************************************************************************
/* 
compile with	
gcc -o rkf45ad rkf45ad.c  -L/usr/local/dislin -ldislnc -lm -lX11
*/

#include "/usr/local/dislin/dislin.h"
#include <stdio.h>
#include <math.h>

#define max 2200

//define r.h.s
double funct_x(double t, double x) {
  return -(sin(pow(t,3))+3*pow(t,3)*cos(pow(t,3)))*x;
}

    double max_steps = 3.0;
    double epstol =  2 * 1E-08;    /* Error control tolerance        */
    int i, j;                    /* Loop counter                   */
    int N=1000;                 /* Number of tentative steps      */
    float  X[max];               /* X array   */
    float T[max];                /* time array */
    double k1, k2, k3, k4, k5, k6;  /* Function values for RKF 4(5)    */
    double s, h;                   /* time step */
    double hmin, hmax;         /* Minimum and maximum step size  */
    double eps, hopt;          /* Error and step size scalar     */
    
    //initial conditions
    double t = 0.0;
    double x = 1.0;

main()
{

 /* Compute the step size */
    h=max_steps/N;
    hmin = h / 32.0;
    hmax = h * 32.0;      
double hin=h;
 printf("Initial time step is h = %lf\n",  h);

  /* Initialize the variable */

    T[0] = 0.;
    X[0] = 1.;
    j    = 0;

 while(T[j] < max_steps)
    {
      if( (T[j] + h) > max_steps ) h = max_steps - T[j];   /* Calculation of the last step */
      t  = T[j];
      x  = X[j];

  /* Compute the function values */

      k1 = h * funct_x(t,x);

      k2 = h * funct_x(t + 0.25 * h, x + 0.25 * k1);

      k3 = h * funct_x(t + 3.0*h/8.0, x + 3.0*k1/32.0 + 9.0*k2/32.0);

      k4 = h * funct_x(t + 12.0*h/13.0, x + 1932.0*k1/2197.0 -
                        7200.0*k2/2197.0 + 7296.0*k3/2197.0);

      k5 = h * funct_x(t + h, x + 439.0*k1/216.0 - 8.0*k2
                         + 3680.0*k3/513.0  - 845.0*k4/4104.0);

      k6 = h * funct_x(t + 0.5*h, x - 8.0*k1/27.0 + 2.0*k2 -
                        3544.0*k3/2565.0 + 1859.0*k4/4104.0 - 11.0*k5/40.0);

      /* Claculate the error eps= |x_(n+1) - x*_(n+1)| */      
      
      eps = fabs( k1/360.0 - 128.0*k3/4275.0 - 2197.0*k4/75240.0
                 + k5/50.0 + 2.0*k6/55.0 );

     if (eps < epstol)
     { X[j+1] = x + 25.0*k1/216.0 + 1408.0*k3/2565.0 + 
                    2197.0*k4/4104.0 - k5/5.0;

       T[j+1] = t + h;
       j++;
    } 

  s=0.84*pow((h*epstol/eps),0.2);

  if (s <= 0.75)
  h /= 2.0;

  else if (s>=2)
  h *= 2.0;

  else
  h=s*h;

  if (h<hmin)
     h=hmin;
  if (h>hmax)
     h=hmax; 

    }
 /* Output */

    for ( i = 0; i <= j; i++)
    {
       printf("i = %d, T[i] = %lf, X[i] = %lf\n", i, T[i], X[i]);
    }
int nn=i-1;
printf("E= %E\n", fabs(exp(-max_steps*sin(pow(max_steps,3)))-X[nn]));
printf("h_in=%lf, h_end= %lf\n", hin, h);

  winsiz (600, 600);
  page (2600, 2600);
  sclmod ("full");
  scrmod ("revers");
  metafl ("xwin");

  disini();
  pagera();
  texmod ("on");
axspos(450,1800);
  axslen(1800,1200);

  name("t","x");
  name("x(t)","y");

  //titlin("$\\exp(-t \\sin(t^3))$",3);

  graf(0.,max_steps,0.,1.0,0.0,18.0,0.,1.0);
  height(50);
  title();

  color("red");
  curve(T,X,nn);
  endgrf();
  disfin();
}