#include<stdio.h>
#include<math.h>
#define size 5 

float  f(float x, float y);
float  Exact(float x);
void   RungeKutta(float x0, float h, float y3[size]);

int main()
{
   int    i;
   float  y1[size], y2[size], y3[size];
   float  h, x0, y0;

   printf("Please enter the step size: ");
    scanf("%f", &h);
   printf("Please enter the initial value x0: ");
    scanf("%f", &x0);
   printf("Please enter the initial value y0: ");
    scanf("%f", &y0);
    
   //Call RungeKutta Method
    y3[0]= y0;
    RungeKutta(x0, h, y3);
     
    //The table
/* y1 -> Approximate values using Euler's Method
   y2 -> Approximate values using Improved Euler Method
   y3 -> Approximate values using Runge Kutta Method
   y  -> Exact solution of the differential equation
   */    
  printf("x \t\t RungeKutta \t Exact \n");
  for(i=0; i < size; i++)
     printf("%f \t %f \t %f \n", x0+i*h, y3[i], Exact(x0+i*h));

    return 0;
}


/* Right Hand Side of the differential Equation 
   y'= -2xy -> f(x,y) = -2xy                    */
float f(float x,float y)
{
  return -2.0*x*y;
}


/* Exact Solution of the given differential equation */
float Exact(float x)
{
  return 2.0*exp(-x*x);
}


/* The Runge-Kutta Method y(x0)=y0 -> 
   k1 = f(x_n,y_n)
   k2 = f(x_n + h/2, y_n +h*k1/2)
   k3 = f(x_n + h/2, y_n +h*k2/2)
   k4 = f(x_(n+1) , y_n +h*k3)

   y_(n+1) = y_n + h*( k1 + 2*k2 + 2*k3 +k4 )/6 */
void RungeKutta(float x0, float h, float y3[size])
{
  int      i=0;
  float   k1, k2, k3, k4, x1;

  while(i < size)
  {

    k1 = f(x0,y3[i]);
    k2 = f(x0+h/2.0, y3[i] + h*k1/2.0);
    k3 = f(x0+h/2.0, y3[i] + h*k2/2.0);
    x1 = x0+h;
    k4 = f(x1, y3[i] + h*k3);

    y3[i+1] = y3[i] + h*(k1 + 2*k2 + 2*k3 + k4)/6.0;
    x0 = x1;
    i++;
  }
}