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

float f1(float x, float y1, float y2)
   {
   return(y2);
   }

float f2(float x, float y1, float y2)
   {
   return( - sin(y1));
   }

int main(void)
{
long int i, j, n;

float h, x0, y10, y20, x, xfim, 
      d1 = 0.0, d2 =0.0, /* erro local estimado. */
      k11, k21, k12, k22, 
      y1, y2, y1p, y2p;

float vf1[2], vf2[2];

x0  = 0.0;
y10 = 0.1;
y20 = 0.0;

h = 0.001;

xfim = 15.0;

n = (int)((xfim - x0)/h);

/* Valores iniciais */
x     = x0;
y1    = y10;
y2    = y20;

vf1[0] = f1(x0, y10, y20);
vf2[0] = f2(x0, y10, y20);

fprintf(stdout, "%e %e %e \n", x, y1, y2, d1, d2);

for (i = 1; i < 2; i++)
  {
   k11    = f1(x, y1, y2);
   k12    = f2(x, y1, y2);
   k21    = f1(x + h, y1 + h * k11, y2 + h * k12);
   k22    = f2(x + h, y1 + h * k11, y2 + h * k12);

   y1     = y1 + (h/2.0) * (k11 + k21);
   y2     = y2 + (h/2.0) * (k12 + k22);

   x      = x + h;

   vf1[i]  = f1(x, y1, y2);
   vf2[i]  = f2(x, y1, y2);

/* ATENCAO!!!!! */
/* Os valores de d1 e d2 não tem significado neste ponto! */
   fprintf(stdout, "%e %e %e %e %e\n", x, y1, y2, d1, d2);
   }

for (i = 2; i < n; i++)
   {
   y1p = y1 + (h/2.0) * (3.0 * vf1[1] - vf1[0]);
   y2p = y2 + (h/2.0) * (3.0 * vf2[1] - vf2[0]);

   x = x + h;

   y1 = y1 + (h/2.0) * (f1(x, y1p, y2p) + vf1[1]);
   y2 = y2 + (h/2.0) * (f2(x, y1p, y2p) + vf2[1]);

   vf1[0] = vf1[1];
   vf2[0] = vf2[1];

   vf1[1] = f1(x, y1, y2);
   vf2[1] = f2(x, y1, y2);

   d1 = -(1.0/6.0) * (y1p - y1);
   d2 = -(1.0/6.0) * (y2p - y2);

   fprintf(stdout, "%e %e %e %e %e\n", x, y1, y2, d1, d2);

   }
}