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

#if !defined(M_E)
#define M_E 2.7182818284590452354
#endif

#define LL 0.0
#define UL 2.0
#define N 200
#define D 5
#define INIT_Y 1.0 /* euler */

double dydx(double x, double y) { return x * y; }

int main() {
  double x, y;
  double y1, y2, h; /* euler 0.01 */
  double y3, y4, h2; /* euler 0.05 */
  double y5, y6, k1, k2, k3, k4; /* Runge-Kutta */
  int i;
  /* euler 0.01 */
  y1 = INIT_Y;
  h =  (UL - LL) / N;
  /* euler 0.05 */
  y3 = INIT_Y;
  h2 = (UL - LL) / N * D;
  /* runge-kutta */
  y5 = INIT_Y;

  for (i = 0; i <= N; i++) {
    x = LL + (UL - LL) / N * i;
    /* true resolution */
    y = exp(0.5 * x * x);
    /* euler 0.01*/
    y2 = y + h * dydx(x, y1);
    /* euler 0.05 */
    if (i % 5 == 0) {
      y4 = y + h2 * dydx(x, y3);
    }
    /* runge-kutta*/
    k1 = dydx(x, y);
    k2 = dydx(x + h / 2.0, y5 + h * k1 / 2.0);
    k3 = dydx(x + h / 2.0, y5 + h * k2 / 2.0);
    k4 = dydx(x + h, y5 + h * k3);
    y6 = y5 + (1.0 / 6.0) * h * (k1 + 2 * k2 + 2 * k3 + k4);
    printf("%.2f, %.10f, %.10f, ", x, y, y1);
    if (i % 5 == 0)
      printf("%.10f, ", y3);
    else
      printf("            , ");
    printf("%.10f\n", y5);
    y1 = y2;
    y3 = y4;
    y5 = y6;
  }
  return 0;
}

  