/*  Approximate the derivative f''(x) using a fourth-order approximation
    of the form R(h) = c1 S(h) + c2 S(2h) where

               35f(x)-104f(x+h/2)+114f(x+h)-56f(x+3h/2)+11f(x+2h)
        S(h) = --------------------------------------------------
                                     3h^2
         

    and constants c1 and c2 are given by Richardson extrapolation. */

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

double f(double x){
    return 1/(x*(x-1)+3);
}
double ddf(double x){
    double t1=2*x-1;
    double t2=x*(x-1)+3;
    double t3=t2*t2;
    return 2*(t1*t1/t2-1)/t3;
}
double S(double x,double h){
    double t=35*f(x)-104*f(x+h/2)+114*f(x+h)-56*f(x+3*h/2)+11*f(x+2*h);
    return t/(3*h*h);
}
double R(double x,double h){
    /* Please fill in this function */
}

#define N 8 

double y[N]= {
  1.50295222758983e-02, 1.09662848519904e-02,
  1.95663737642718e-03, 1.48195086357439e-04,
  8.54554946286279e-06, 4.73880335555288e-07,
  2.70769011656569e-08, 1.58625496071707e-09
};

int main(){
    printf("This is problem 2 on final part 2.\n");
    printf("%21s %21s %21s\n","h","R_h","Error");
    double eold=1,x=1,ddfx=ddf(x),h=1,flag=0;
    for(int n=0;n<N;n++){
        double t=R(x,h),e=fabs(t-ddfx);
        if(fabs(e-y[n])>y[n]*1e-6) flag=1;
        printf("%21.14e %21.14e %21.14e\n",h,t,e);
        h/=2;
    }
    if(!flag){
        printf("Congratulations!  Your answer is correct!\n");
    } else {
        printf("Please try again.  Your answer is incorrect.\n");
    }
    return 0;
}