#include <stdio.h>
#include <math.h>

double f(double x,double y){
    return x*cos(y*x);
}

static double h=0.1;
double rk2step(double x,double y){
    double k1=h*f(x,y);
    return y+(k1+h*f(x+h,y+k1))/2;
}
// Runge-Kutta Order Four from page 288
double rk4step(double t,double w){
    double k1=h*f(t,w);
    double k2=h*f(t+h/2,w+k1/2);
    double k3=h*f(t+h/2,w+k2/2);
    double k4=h*f(t+h,w+k3);
    return w+(k1+2*(k2+k3)+k4)/6;
}

int main(){
    FILE *output=fopen("adapt.dat","w");
    double x0=0,xn=10;
    double y0=2,tol=1e-6;
    double x=x0,y=y0,esofar=0;
    fprintf(output,"%g %g\n",x,y);
    for(int fstop=0;fstop!=1;){
        if(x+h>xn) {
            h=xn-x;
            fstop=1;
        }
        double t1=rk2step(x,y);
        double t2=rk4step(x,y);
        double epsilon=h/(xn-x)*(tol-esofar);
        double steperr=fabs(t1-t2);
        if(steperr>epsilon){
            double alpha=pow(0.9*epsilon/steperr,1.0/3);
            h*=alpha;
            printf("at %g step size is %g\n",x,h);
            fstop=0;
        } else {
            y=t2; x+=h;
            esofar+=steperr;
            if(steperr<epsilon/3){
                h*=pow(2.0,1.0/3);
                printf("step bigger now!\n");
            }
        }
        fprintf(output,"%g %g\n",x,y);
    }
    fclose(output);
    return 0;
}