#include <stdio.h>

double sigma = 10.0;
double rho = 28.0;
double beta = 8./3.;

double fx(const double x, const double y, const double z) {
	return sigma*(y-x);
}

double fy(const double x, const double y, const double z) {
	return x*(rho-z)-y;
}

double fz(const double x, const double y, const double z) {
	return x*y-beta*z;
}

int main() {
	double tmax = 50.0;
	double dt = 0.001;
	int N = tmax/dt;

	double x,y,z,t;

	//initial conditions, i.e. x(0),y(0),z(0)
	x = 0.0;
	y = 2.0;
	z = 20.0;

	// open file
	// this file would contain 3 columns
	// t, x, y, z
	FILE *fp = fopen("butterfly.txt","w");

	// required for RK4
	// see the iteration rule
	double k1,k2,k3,k4;
	double m1,m2,m3,m4;
	double n1,n2,n3,n4;

	for(int i=0; i<N-1; i++) {

		t = i*dt;

		//write to file
		fprintf(fp,"%f %f %f %f\n",t,x,y,z);


		k1 = dt*fx(x,y,z);
		m1 = dt*fy(x,y,z);
		n1 = dt*fz(x,y,z);

		k2 = dt*fx(x+k1/2, y+m1/2, z+n1/2);
		m2 = dt*fy(x+k1/2, y+m1/2, z+n1/2);
		n2 = dt*fz(x+k1/2, y+m1/2, z+n1/2);

		k3 = dt*fx(x+k2/2, y+m2/2, z+n1/2);
		m3 = dt*fy(x+k2/2, y+m2/2, z+n2/2);
		n3 = dt*fz(x+k2/2, y+m2/2, z+n2/2);

		k4 = dt*fx(x+k3, y+m3, z+n1);
		m4 = dt*fy(x+k3, y+m3, z+n3);
		n4 = dt*fz(x+k3, y+m3, z+n3);

		x = x + (k1+(k2*2)+(k3*2)+k4)/6.0;
		y = y + (m1+(m2*2)+(m3*2)+m4)/6.0;
		z = z + (n1+(n2*2)+(n3*2)+n4)/6.0;

	}


	return 0;
}