ホイン法(heun method)

/**********************************************************************
y' = x+ y, 0<=x<=1,y(0)=1の常微分方程式の解をもとめる
 **********************************************************************/
#include <stdio.h>
#define N 20

double f1(double x,double y);
void heun(double x,double y,double a,double b, int n, double(*f)(double,double));

int main(void)
{
  heun(0.0, 1.0, 0.0, 1.0, N, f1 );
  return 0;
}
void heun(double x,double y,double a,double b,int n,double (*f)(double,double))
{
  int i;
  double h,k1,k2;
  h = (b-a)/n;
  x = a;
  for(i =0; i<=n-1; i++){
    k1 = f(x,y);
    k2 = f(x + h, y + h*k1 );
    y += h / 2 * ( k1 + k2);
    x += h;
    printf("x=%lf y=%lf\n",x,y);
  }
}

double f1(double x,double y)
{
  return x+y;
}

/**********************************************************************
$ ./a.out
x=0.050000 y=1.052500
x=0.100000 y=1.110253
x=0.150000 y=1.173529
x=0.200000 y=1.242609
x=0.250000 y=1.317793
x=0.300000 y=1.399393
x=0.350000 y=1.487736
x=0.400000 y=1.583170
x=0.450000 y=1.686058
x=0.500000 y=1.796781
x=0.550000 y=1.915741
x=0.600000 y=2.043360
x=0.650000 y=2.180082
x=0.700000 y=2.326374
x=0.750000 y=2.482726
x=0.800000 y=2.649653
x=0.850000 y=2.827698
x=0.900000 y=3.017430
x=0.950000 y=3.219448
x=1.000000 y=3.434382
 **********************************************************************/