#include <stdio.h>
#include <math.h>
/* adm.c : 1999 3 28 by T.Oguni */
#define N 200
double f(double t, double y) {return 2.*t*y; }
adm(tvals,yvals,t0,tf,y0,h)
double tvals[N+1], yvals[N+1], t0, tf, y0, h;
{
  double k[4],fval[5],ys[4],y,t,s,y1,t1,yc;
  double b[4]={0.1666666667, 0.333333333, 0.3333333, 0.166666667};
  double d[4]={0., 0.5, 0.5, 1.};
  double c[4][4]={0.,0.,0.,0., 0.5,0.,0.,0., 0.,0.5,0.,0., 0.,0.,1.,0.};
  int order=4, i, j, m, n;
  y=y0; t=t0; yvals[0]=y0; tvals[0]=t0;
  fval[0]=f(t,y); ys[0]=y0; s=(tf-t0)/h; n=(int)s+1;
  for (j=1; j<4; j++){
    k[0]=h*f(t,y);
    for (i=1; i<order; i++){
      s=0.;
      for (m=0; m<i; m++){
        s+=c[i][m]*k[m]; }
      k[i]=h*f(t+h*d[i], y+s);
    }
    s=0.;
    for (m=0; m<order; m++){
      s+=b[m]*k[m]; }
    y1=y+s; ys[j]=y1; t1=t+h;
    fval[j]=f(t1,y1); tvals[j]=t1; yvals[j]=y1;
    t=t1; y=y1;
  }
  for (i=4; i<n; i++){
    y1=ys[3]+h*(55.*fval[3]-59.*fval[2]+37.*fval[1]-9.*fval[0])/24.;
    t1=t+h; fval[4]=f(t1,y1);
    yc=ys[3]+h*(9.*fval[4]+19.*fval[3]-5.*fval[2]+fval[1])/24.;
    fval[4]=f(t1,yc);
    for (m=0; m<4; m++){ fval[m]=fval[m+1]; }
    ys[3]=yc; tvals[i]=t1; yvals[i]=yc;
    t=t1; y=y1;
  }
  return;
}
/* main */
main(){
  double tvals[N+1],yvals[N+1],t0,tf,y0,h,s;
  int i, n;
  t0=0.; tf=2.; h=0.25; y0=2.; s=(tf-t0)/h; n=(int)s+1;
  adm(tvals,yvals,t0,tf,y0,h);
  for (i=0; i<n; i++){
    printf("%f %f\n", tvals[i], yvals[i]); }
}