3. 3章の答え・２

#include <stdio.h>
#include <math.h>
#define N 100
#define Xs 0
#define Xe 3.0
double f(double x, double y);
double analytical(double x);
int main(void) {
  double h=(Xe-Xs)/N;
  double x=0.0;
  double y=4.0;
  double k1,k2,k3,k4;
  int i;
  for (i = 0; i < N; i++) {
    k1=h*f(x, y);
    k2=h*f(x+h/2.0, y+k1/2.0);
    k3=h*f(x+h/2.0, y+k2/2.0);
    k4=h*f(x+h, y+k3);
    y=y+(k1+2.0*k2+2.0*k3+k4)/6.0;
    x=x+h;
    printf("%4.4lf\t %4.4lf\t  %4.4lf\n",x,y,analytical(x));   
  }
}
double f(double x, double y) {
        return x * y;
}
double analytical(double x) {
        return 4.0*exp(x*x / 2.0);
}