/* RK4.c   -- solution of IVP by the Runge Kutta Method
 this program provides an approximate (numerical) solution to
 the well-posed initial value problem y' = f(t,y), a <= t <= b, (1)
 y(a) = alpha at (n+1) equally-spaced abscissa points in the interval
 [a,b].  note that y is taken as the dependent function of the
 independent variable t.

 The Runge Kutta Method of Order 4 constructs an approximation w[i]
 to the true solution
 y(t[i]) at each mesh point in the interval [a,b] for each
 i = 1,2, ... , n by the recursive scheme

 w[0] = alpha
 w[i+1] = w[i] + (K1 + 2K2 + 2K3 + K4)/6   */

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

double func(double t,double y);      //function prototype for y ' = f(t,y)

int main(void)
{
   double a,b,alpha,h,t,w,K1,K2,K3,K4;
   int n;
   int i;
   printf("RungeKutta4 Method for solution of y' = f(t,y) on [a,b], y(a)=alpha\n");
   printf("Please enter interval limits [a,b], b > a:\n");
   scanf("%lf %lf",&a,&b);
   while (b<=a)
   {
      printf("b must be greater than a.  Please re-enter a and b: \n");
      scanf("%lf %lf",&a,&b);
   }
   printf("Please enter number of steps n to take in [a,b]:\n");
   scanf("%d",&n);
   while (n<=0)
   {
	  printf("n must be a positive integer. Please re-enter n: \n");
      scanf("%d",&n);
   }

   h = (b-a)/n;
   printf("Please enter initial condition alpha, y[%lf] = alpha:\n",a);
   scanf("%lf",&alpha);
   /*  INITIALIZATION       */
   i = 0;
   t=a;
   w = alpha;
   /* END OF INITIALIZATION */
   printf("\n\nIteration\t\tt\t\t\ty\n");
   printf("%d\t\t\t%6.6lf\t\t%6.6lf\n",i,t,w);
   for (i = 1;i<=n;i++)
   {
	   K1 = h*func(t,w);
	   K2 = h*func(t+h/2,w+K1/2);
	   K3 = h*func(t+h/2,w+K2/2);
	   K4 = h*func(t+h,w+K3);
	   w = w + (K1 + 2*K2 + 2*K3 + K4)/6.;
	   t = a + i*h;
	   printf("%d\t\t\t%6.6lf\t\t%6.6lf\n",i,t,w);
   }
   return 0;
}

double func(double t,double y)
{
	return -y + t + 1;
}