#include <math.h>
#include <stdio.h>


#define F(t)	((t) + 1.0 / (1.0 - (t)))
#define f(t, u)	(1.0 + ((t) -(u)) * ((t) - (u)))


#define t0	((double) 2.0)
#define t1 	((double) 3.0)
#define dt 	((double) 1.0e-8)
#define u0 	((double) 1.0)
#define n	(1 + (int) ((t1 - t0) / dt))


double t[n];
double y[n];
double u[n];


void init()
{
	int i;
	
	t[0] = t0;
	y[0] = F(t0);

	for (i=1; i<n; ++i)
	{
		t[i] = t[i-1] + dt;
		y[i] = F(t[i]);
	}
}


void euler()
{
	int i;

	u[0] = u0;

	for (i=1; i<n; ++i)
		u[i] = u[i-1] + dt * f(t[i-1], u[i-1]);
}


void runge_kutta()
{
	int i;

	u[0] = u0;

	for (i=1; i<n; ++i)
	{
		double k0 = dt * f(t[i-1], u[i-1]);
		double k1 = dt * f(t[i-1] + dt/2.0, u[i-1] + k0/2.0);
		double k2 = dt * f(t[i-1] + dt/2.0, u[i-1] + k1/2.0);
		double k3 = dt * f(t[i], u[i-1]+k2);
		u[i] = u[i-1] + (k0 + 2.0*k1 + 2.0*k2 + k3) / 6.0;
	}
}

void print()
{
	int i;

	printf("i    \tt                     \ty                     \tu                     \te                     \n"
	       "-----\t----------------------\t----------------------\t----------------------\t----------------------\n");

	for (i=0; i<n; ++i)
		printf("%05d\t%+.15e\t%+.15e\t%+.15e\t%+.15e\n", i, t0+i*dt, y[i], u[i], y[i]-u[i]);
}


int main()
{
	init();
	euler();
	print();
	runge_kutta();
	print();
	return 0;
}