#include <stdio.h>
#include <math.h>
/* compile using the command

    gcc rk4.c -lm

More information on Taylor method and the
Trapezoid method can be found in Burden,
Faires and Burden, Numerical Analysis 10th edition
Chapter 5.
*/

static double h=1.0/1024;
double exacty(double x){
    return x*x*cos(x);
}
double f(double x,double y){
    return 2/x*y-x*x*sin(x);
}
double ginv(double x,double g){
    return (g-h*x*x/2*sin(x))/(1-h/x);
}
double fx(double x,double y){
    return -2*y/(x*x)-2*x*sin(x)-x*x*cos(x);
}
double fy(double x,double y){
    return 2/x;
}
double trapezoidstep(double x,double y){
    return ginv(x+h,y+h/2*f(x,y));
}
double taylorstep(double x,double y){
    double t=f(x,y);
    return
    y+h*(t+h/2*(fx(x,y)+fy(x,y)*t));
}
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;
}
typedef double (*odesolver)(double,double);
int main(){
    double x,Y[4];
    odesolver step[4]={taylorstep,trapezoidstep,
        rk2step,rk4step};
    char *names[4]={"taylor","trapezoid","rk2","rk4"};
    for(int n=2;n<=100000;n*=2){
        h=1.0/n;
        for(int k=0;k<4;k++) Y[k]=cos(1);
        int j;
        for(j=0;j<n;j++){
            x=1+j*h;
            for(int k=0;k<4;k++) Y[k]=step[k](x,Y[k]);
        }
        x=1+j*h;
        double ye=exacty(x);
        printf("%9d ",n);
        for(int k=0;k<4;k++)
            printf("%16.9e ",fabs(Y[k]-ye));
        printf("\n");
    }
    return 0;
}