------------------------------------------------------------------

#include<stdio.h> 

#include<math.h> 

double f1(double t, double y);//문제 식 정의

double f2(double t, double y);

double f3(double t, double y, double z );

void main()

{

double y1[200];//y1를 저장할 공간 

double k[100][100];//k값을 저장할 공간 

y1[0] = 0;//y1의 초기치

double y1p, y1c, y1pm;// y1의 예측자, 수정자, 보정예측자 

double t = 0;

int i, j;

double h = 0.1;

for (i = 0; i<3; i++)

{

k[1][1] = f1(t, y1[i]);

k[1][2] = f1(t + h / 2, y1[i] + h / 2.0 * k[1][1]);

k[1][3] = f1(t + h / 2, y1[i] + h / 2.0 * k[1][2]);

k[1][4] = f1(t + h, y1[i] + h * k[1][3]);//y1의 4계 Runge-Kutta식 

y1[i + 1] = y1[i] + h / 6.0 * (k[1][1] + 2 * (k[1][2] + k[1][3]) + k[1][4]);

}

for (j = 3; j<100; j++)

{

t = 0.1 + j / 10.0;

y1p = y1[j] + h / 24.0*(55 * f1(t - 0.1, y1[j]) - 59 * f1(t - 0.2, y1[j - 1]) + 37 * f1(t - 0.3, y1[j - 2]) - 9 * f1(t - 0.4, y1[j - 3]));

y1c = y1[j] + h / 24.0*(9 * f1(t, y1p) + 19 * f1(t - 0.1, y1[j]) - 5 * f1(t - 0.2, y1[j - 1]) + f1(t - 0.3, y1[j - 2]));

//y1의 예측자와 수정자를 구하는 식 

y1pm = y1p + 251.0 / (19.0 + 251.0)*(y1c - y1p);

//y1의 보정 예측자를 구하는 식 

y1c = y1[j] + h / 24.0*(9 * f1(t, y1pm) + 19 * f1(t - 0.1, y1[j]) - 5 * f1(t - 0.2, y1[j - 1]) + f1(t - 0.3, y1[j - 2]));

y1[j + 1] = y1c - 19.0 / (19.0 + 251.0) * (y1c - y1p);

//y1보정 수정자를 구하는 식 

if (t > 10)break;

printf(" t=%lf, y1=%lf \n", t, y1[j + 1]);

}

}

double f1(double t, double y, double z)

{

double f;

f = y ;

return f;

}

double f2(double t, double y, double z)

{

double f;

f = z ;

return f;

}

double f3(double t, double y, double z)

{

double f;

f = 2*t+3*y-t*z;

return f;

}