// Runge-Kutta order 4

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

#define F(t,y)      (-(y+1.0)*(y+3.0))
#define ExactF(t)   (-3.0+2.0/(1.0+exp(-2*(t))))   
#define t0        0.0   // tÀÇ ÁÂÃø °æ°èÄ¡
#define te        2.0   // tÀÇ ¿ìÃø °æ°èÄ¡  
#define y0       -2.0   // yÀÇ ÃʱâÄ¡
#define num      20.0   // µîºÐ¼ö

void main()
{
   double k1, k2, k3, k4, t, w[100];
   double h;
   int i;
   char ExFSel;
  
   t = t0; w[0] = y0;  // ÃʱâÄ¡ ´ëÀÔ
   h = (te-t0) / num;

   for(i=0; i<num; i++) {
      k1 = h * F(t, w[i]);
      k2 = h * F(t+h/2.0, w[i] + k1/2.0);
      k3 = h * F(t+h/2.0, w[i] + k2/2.0);
      k4 = h * F(t+h, w[i] + k3);

      t+=h;

      w[i+1] = w[i] + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
   }

   // Á¤È®ÇÑ ÇÔ¼ö°¡ ÀÖ´ÂÁö ¿©ºÎ
   printf("Does Exact Function Exists?");
   scanf("%c", &ExFSel);

   if(ExFSel=='y' | ExFSel=='Y') {
      printf("ti         Exact             wi            Error\n");        
      for(i=0; i<=num; i++)
         printf("%4.2f  %12.8f     %12.8f     %12.8f\n",
              t0+h*i, ExactF(t0+h*i), w[i], fabs(ExactF(t0+h*i)-w[i]));
   }          
   else {
      printf("ti             wi\n");
      for(i=0; i<=num; i++)
         printf("%4.2f       %12.8f\n", t0+h*i, w[i]);
   }
   
}