#include <iostream>
#include <cmath>
double f(double, double, double); // Интерфейс функции f, это функция при x'
double g(double, double, double); // Интерфейс функции g, это при y'
double distance(double, double, double*);// Уже не нужен
using namespace std;// Это чтоб проще было пользоваться вводом-выводом
int main() {
  int k = 0, c1, c2;//k - переменная-счетчик количества точек в графике, используется для ограничения, c1 и c2 вводились для построения множества графиков, сейчас вместо них i и j
  double t = 0, X, Y, h = 0.01;// t - время, X и Y - значения функций, h - шаг моделирования
  double k1, k2, k4, k3, q1, q2, q4, q3;//служебные переменные
  double X0 = -1.5 * M_PI, Y0 = 0;//С помощью X0 и Y0 задаются координаты особой точки для удобства
  std::cout.precision(10);//Задается точность вывода чисел в потоке stdout
  
  for(double j = -100; j <= 100; j+=1)
    for(double i = -10; i <= 10; i+=1) {//Циклы, двигающие начальную точку. Можно их убрать, тогда будет строиться одна траектория
      X = j*0.1;
      Y = i*0.1;//Задание начальной точки моделирования
      t=0;
      k=0;
      cout << "#X=" << X << " Y=" << Y << "\n";//Вывод первичного значения для понимания, откуда траетория начиналась, # позволяет gnuplot'у игнорировать эту строку
      for(; (fabs(X) <= 10 && fabs(Y) <= 10)  && k++<1e4; t += h){//Цикл построения траектории. Условие ограничивает область [-10:10] по X и [-10:10] по Y и не более 10000 точек на траекторию
	cout<< t << ' ' << X << ' ' << Y << "\n";// Вывод координат точки
	k1 = h * g(t, X, Y);
	q1 = h * f(t, X, Y);
	k2 = h * g(t + h/2.0, X + q1/2.0, Y + k1/2.0);
	q2 = h * f(t + h/2.0, X + q1/2.0, Y + k1/2.0);
	k3 = h * g(t + h/2.0, X + q2/2.0, Y + k2/2.0);
	q3 = h * f(t + h/2.0, X + q2/2.0, Y + k2/2.0);
	k4 = h * g(t + h, X + q3, Y + k3);
	q4 = h * f(t + h, X + q3, Y + k3);
	Y = Y + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
	X = X + (q1 + 2.0*q2 + 2.0*q3 + q4)/6.0;
      }
      cout<<"\n";//Вывод пстой строки для удобства построения графиков
    }
}
double f(double t, double x, double y){
  return y;
}
double g(double t, double x, double y){
  return 2.01 * y - cos(x);
}