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

    gcc trapezoid.c -lm

plot output using gnuplot and the commands

    set logscale x
    set logscale y
    plot "trapezoid.dat" ti "taylor", \
    "" using 1:3 ti "trapezoid",1/x**2

More information on Taylor method and the
Trapezoid method can be found in Burden,
Faires and Burden, Numerical Analysis
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;
}
int main(){
    double x,ytaylor,ytrapezoid,yrk2;
    printf("#%8s %16s %16s %16s\n","n","error-taylor",
        "error-trapezoid","error-rk2");
    for(int n=64;n<=100000;n*=2){
        h=1.0/n;
        ytaylor=cos(1); ytrapezoid=ytaylor;
        yrk2=ytaylor;
        int j;
        for(j=0;j<n;j++){
            x=1+j*h;
            ytaylor=taylorstep(x,ytaylor);
            ytrapezoid=trapezoidstep(x,ytrapezoid);
            yrk2=rk2step(x,yrk2);
        }
        x=1+j*h;
        double ye=exacty(x);
        printf("%9d %16.9e %16.9e %16.9e\n",n,
            fabs(ytaylor-ye),fabs(ytrapezoid-ye),
            fabs(yrk2-ye));
    }
    return 0;
}