/**
 * numerical solution of y'=f(x,y), y(x0)=y0
 **/

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

/**
 * the integration function:  f(x,y) = y
 **/
double function(double x, double y) {
  return y;
}

/**
 * a single step solver:  Euler method
 * y(j+1) = y(j) + h f(xj,yj)
 **/
double step(double y0, double x0, double h) {
  return y0 + h * function(x0, y0);
}

/**
 * exact solution
 **/
double exact(double x) {
  return exp(x);
}

/**
 * print out one result row.
 * compare with exact (less approximate) value.
 **/
void dumprow(double x, double y, double yexact) {
  printf("%10lf %10lf %10lf %10lf\n", x, y, yexact, yexact-y);
}

/**
 * main function
 **/
int main(int argc, char *argv[]) {

  double x,y,xmax,xstep;

  /* take input */
  printf("Input the maximum:  ");
  scanf("%lf", &xmax);
  printf("Input the increment:  ");
  scanf("%lf", &xstep);

  /* initialize variables */
  x = 0.0;
  y = 1.0;
  dumprow(x, y, exact(x));

  while (x < xmax) {

    y = step(y, x, xstep);
    x += xstep;

    /* print result */
    dumprow(x, y, exact(x));
  }

  return 0;
}