#include <iostream>
#include <iomanip>

#include <math.h>

typedef double (*ODE_RHS) (double, double);

double RightHandSide(double t, double x)
{
   return (t*x-x*x)/(t*t);
}

/*
  A function implementing the second order Runge-Kutta method for
  solving the ordinary differential equation
                       x'=f(t,x), x(t0)=x0
  in the interval t0 -> tn, with step-size h=(tn-t0)/n.
  The solution is returned in the vector x (of length n+1 (!)),
  where x[i] = x(t0+i*h) (and consequently x[0]=t0, x[n]=tn).
*/  
double runge_kutta_2(ODE_RHS f, double t0, double tn,
                     int n, double x0, double *x)
{
   double h = (tn-t0)/n;
   x[0] = x0;

   for (int i=1; i<=n; i++)
   {
      double t = t0+(i-1)*h, xt=x[i-1];  // shorthand notations
      double F1 = h*f(t,xt);             // F1 = hf(t,x)
      double F2 = h*f(t+h,xt+F1);        // F2 = hf(t+h,x+F1)
      x[i] = xt + 0.5*(F1+F2);           // x(t+h) = x(t) + 0.5(F1+F2)
   }
}

/****************************************************************/

double solution(double t)
{
   return t/(0.5+log(t));
}

int main()
{
   int N=256;
   double *x = new double[N+1];

   runge_kutta_2(RightHandSide,1.0,3.0, N,2.0,x);
   
   std::cout << std::setiosflags(std::ios::fixed) << std::setprecision(8);
   for(int i=0; i<=N; i++)
   {
      std::cout << "t=" << 1+i*(2.0/N) << ", ~x(t)="
                << x[i] << ", x(t)=" << solution(1+i*(2.0/N)) << std::endl;
   }

   delete x;
}