#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>

using namespace std;

double fun1(double yrab, double yfox)
{
	return 1.2e0*yrab-0.6e0*yrab*yfox;
}
double fun2(double yrab, double yfox)
{
	return -0.8e0*yfox+0.3e0*yrab*yfox;
}

int main(int argc, char** argv)
{
	vector <double> y(2),aux(2);
	vector < vector <double> > ak(2, vector <double> (4));

	double time=0.e0;
	double deltat=0.1;

	int neqs=2;
	y[0]=10;
	y[1]=2;
	cout<<time<<"\t";
	for(int i=0;i<neqs;i++)
		cout<<y[i]<<"\t";
	cout<<endl;
	
	for(int ktimestep=1;ktimestep<500;ktimestep++)
	{
	//Step 1:
		ak[0][0] = deltat*fun1(y[0],y[1]);
		ak[1][0] = deltat*fun2(y[0],y[1]);
		aux[0] = y[0]+0.5e0*ak[0][0];
		aux[1] = y[1]+0.5e0*ak[1][0];
       //Step 2:
		ak[0][1] = deltat*fun1(aux[0],aux[1]);
		ak[1][1] = deltat*fun2(aux[0],aux[1]);
		aux[0] = y[0]+0.5e0*ak[0][1];
		aux[1] = y[1]+0.5e0*ak[1][1];
	//Step 3:
		ak[0][2] = deltat*fun1(aux[0],aux[1]);
		ak[1][2] = deltat*fun2(aux[0],aux[1]);
		aux[0] = y[0]+ak[0][2];
		aux[1] = y[1]+ak[1][2];
	//Step 4:
		ak[0][3] = deltat*fun1(aux[0],aux[1]);
		ak[1][3] = deltat*fun2(aux[0],aux[1]);
	//Synthesis:
		y[0]=y[0]+(ak[0][0]+2.e0*ak[0][1]+2.e0*ak[0][2]+ak[0][3])/6.e0;
		y[1]=y[1]+(ak[1][0]+2.e0*ak[1][1]+2.e0*ak[1][2]+ak[1][3])/6.e0;
		time=time+deltat;

		cout<<time<<"\t";
		for(int i=0;i<neqs;i++)
			cout<<y[i]<<"\t";
		cout<<endl;
	}
	return 0;
}