RK

//IEulerRK.cpp--Improved Euler and Runge-Kutta(4)
//qiu changweifen fangcheng shuzhijie
#include < stdio.h > 
#include < math.h > 
#define FMT "%-15.7g"
typedef double dbl;
//prototypes
dbl fxy(dbl x, dbl y);
dbl f(dbl x);
void RK(dbl x0, dbl y0, dbl b, dbl h);
int main() {
    printf("\nRunge-Kutta:\nh=0.1:\n");
    RK(0.0, 1.0 / 3.0, 1.0, 0.1);
    printf("\nh=0.025:\n");
    RK(0.0, 1.0 / 3.0, 1.0, 0.025);
    printf("\nh=0.01:\n");
    RK(0.0, 1.0 / 3.0, 1.0, 0.01);
    return 0;
}
dbl fxy(dbl x, dbl y) {
    return ( - 50 * y + 50 * x * x + 2 * x);
}
dbl f(dbl x) {
    return (exp( - 50 * x) / 3.0 + x * x);
}
void RK(dbl x0, dbl y0, dbl b, dbl h) {
    dbl k1,
    k2,
    k3,
    k4;
    dbl xn,
    yn;
    printf("xn             yn             y(xn)          \n");
    xn = x0;
    yn = y0;
    while (xn <= b - 1e-7) {
        k1 = fxy(xn, yn);
        k2 = fxy(xn + h / 2.0, yn + k1 * h / 2.0);
        k3 = fxy(xn + h / 2.0, yn + k2 * h / 2.0);
        k4 = fxy(xn + h, yn + k3 * h);
        yn = yn + (k1 + 2 * k2 + 2 * k3 + k4) * h / 6;
        xn += h;
        printf(FMT, xn);
        printf(FMT, yn);
        printf(FMT, f(xn));
        printf("\n");
    }
}