#include <stdio.h>

double f(double x, double y) { return (y * y - x) / (y + x * y * y); };

double rungen_kutta_4(double x0, double y0, double xn, int steps) {
  double h = (xn - x0) / steps;

  while (steps--) {
    double k1 = f(x0, y0) * h;
    double k2 = f(x0 + h / 2, y0 + k1 / 2) * h;
    double k3 = f(x0 + h / 2, y0 + k2 / 2) * h;
    double k4 = f(x0 + h, y0 + k3) * h;

    double k = (k1 + k4 + 2 * (k2 + k3)) / 6;

    x0 += h, y0 += k;
  }

  return y0;
}

int main() {
  double x0 = 0, y0 = 5, xn = 1784;
  int steps = 1000000;

  printf("Value of y at x=%0.4lf is y=%0.4lf", xn,
         rungen_kutta_4(x0, y0, xn, steps));
}