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

double t=0;
double y=1;
double T=3;
double error=0.01;
double f(double t, double y){
    return y-y*sin(t*y);
}

int main(){
    double k1,k2,k3,k4,k5,k6,y1,y2;
    double h=(T-t)/10;
    for(;t<T;){
        if(t+h>T) h=T-t;
        k1=f(t,y);
        k2=f(t+h/4,y+h/4*k1);
        k3=f(t+h*3/8,y+3*h/32*k1+9*h/32*k2);
        k4=f(t+h*12/13,y+1932*h/2197*k1-7200*h/2197*k2+7296*h/2197*k3);
        k5=f(t+h,y+439*h/216*k1-8*h*k2+3680*h/513*k3-845*h/4104*k4);
        k6=f(t+h/2,y-8*h/27*k1
            +2*h*k2-3544*h/2565*k3+1859*h/4104*k4-11*h/40*k5);
        y1=y+h*(16*k1/135+6656*k3/12825+28561*k4/56430-9*k5/50+2*k6/55);
        y2=y+h*(25*k1/216+1408*k3/2565+2197*k4/4104-k5/5);
        double e=fabs(y1-y2);
        if(e<error*h){
            y=y1;
            t=t+h;
            printf("%g %g\n",t,y);
            if(128*e<error*h) {
                printf("# step too small--doubling next one!\n");
                h=2*h;
            }
        } else {
            h=h/2;
            printf("# step too big--halve and try again!\n");
        }
    }
    return 0;
}