#include <bits/stdc++.h>
using namespace std;

int main() {
  auto dy1 = [&] (double y1, double y2) { return 500 * y2 - 501 * y1; };
  auto dy2 = [&] (double y1, double y2) { return 500 * y1 - 501 * y2; };

  double h = 0;
  cout << "Write your precesion" << endl;
  cin >> h;


  double y1_0 = 5;
  double y2_0 = -10;

  cout << "Write initial conditions" << endl;
  cin >> y1_0;
  cin >> y2_0;

  vector<double> y1s, y2s;
  y1s.push_back(y1_0);
  y2s.push_back(y2_0);

  vector<double> coefficients = { 55. / 24, -59. / 24, 37. / 24, -09. / 24 };
  int cnt = 0;
  // Eulers method
  cout << "Getting start 4 points unsing Eulers method" << endl;
  for (int i = 1; i < 4; i++) {
    double nextY1 = y1s.back() + h * dy1(y1s.back(), y2s.back());
    double nextY2 = y2s.back() + h * dy2(y1s.back(), y2s.back());
    y1s.push_back(nextY1);
    y2s.push_back(nextY2);
    cnt++;
    cout << fixed << setprecision(10) << "Iterations: " << cnt << " y1: " << y1s.back() << " y2: " << y2s.back() << endl;
  }

  cout << "Write number of iterations" << endl;
  int n;
  cin >> n;
  auto getNext = [&] (int i, const function<double(double, double)>& f, double y) {
    double dlt = 0;
    for (int j = 1; j <= 4; j++) {
      dlt += coefficients[j - 1] * f(y1s[i - j], y2s[i - j]);
    }
    return y + h * dlt;
  };
  for (; cnt <= n; cnt++) {
    double nextY1 = getNext(cnt, dy1, y1s.back());
    double nextY2 = getNext(cnt, dy2, y2s.back());
    y1s.push_back(nextY1);
    y2s.push_back(nextY2);
    cout << fixed << setprecision(10) << "Iterations: " << cnt << " y1: " << y1s.back() << " y2: " << y2s.back() << endl;
  }

}