#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
double *dvector(long i, long j);
void free_dvector(double *a);
void FUNC(double x, double *y, double *f);
void rk4m(double *y, double *f, int N, double a, double b, double step, void (*F)(double, double *, double *));
 
int main(void) {
    int N = 2;
    double *y, *f, a = 0.0, b = 1.0, step = 0.01;
 
    y = dvector(0, N - 1);
    f = dvector(0, N - 1);
 
    y[0] = 0.1; 
    y[1] = 0.0;
 
    printf("x\ty\n"); // ヘッダーを出力
    rk4m(y, f, N, a, b, step, FUNC); // ファイル操作をせず標準出力
 
    free_dvector(y);
    free_dvector(f);
 
    return 0;
}
 
void FUNC(double x, double *y, double *f) {
    f[0] = y[1]; // dy1/dx = y2
    f[1] = cos(x) - 4 * y[1] - 3 * y[0]; // dy2/dx = cos(x) - 4y2 - 3y1
}
 
void rk4m(double *y, double *f, int N, double a, double b, double step, void (*F)(double, double *, double *)) {
    double *k1, *k2, *k3, *k4, *tmp, x, h;
    int j;
 
    k1 = dvector(0, N - 1);
    k2 = dvector(0, N - 1);
    k3 = dvector(0, N - 1);
    k4 = dvector(0, N - 1);
    tmp = dvector(0, N - 1);
 
    h = step;
    x = a;
 
    while (x <= b) {
        printf("%.2lf\t%.10lf\n", x, y[0]); // x と y[0] を標準出力
 
        F(x, y, f);
        for (j = 0; j < N; j++) k1[j] = f[j];
 
        for (j = 0; j < N; j++) tmp[j] = y[j] + h * k1[j] / 2.0;
        F(x + h / 2.0, tmp, f);
        for (j = 0; j < N; j++) k2[j] = f[j];
 
        for (j = 0; j < N; j++) tmp[j] = y[j] + h * k2[j] / 2.0;
        F(x + h / 2.0, tmp, f);
        for (j = 0; j < N; j++) k3[j] = f[j];
 
        for (j = 0; j < N; j++) tmp[j] = y[j] + h * k3[j];
        F(x + h, tmp, f);
        for (j = 0; j < N; j++) k4[j] = f[j];
 
        for (j = 0; j < N; j++)
            y[j] = y[j] + h * (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6.0;
 
        x += h; // x を更新
    }
 
    free_dvector(k1);
    free_dvector(k2);
    free_dvector(k3);
    free_dvector(k4);
    free_dvector(tmp);
}
 
double *dvector(long i, long j) {
    double *a;
    if ((a = malloc((j - i + 1) * sizeof(double))) == NULL) {
        printf("メモリが確保できません\n");
        exit(1);
    }
    return a;
}
 
void free_dvector(double *a) {
    free(a);
}
 

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KI2luY2x1ZGUgPG1hdGguaD4KCmRvdWJsZSAqZHZlY3Rvcihsb25nIGksIGxvbmcgaik7CnZvaWQgZnJlZV9kdmVjdG9yKGRvdWJsZSAqYSk7CnZvaWQgRlVOQyhkb3VibGUgeCwgZG91YmxlICp5LCBkb3VibGUgKmYpOwp2b2lkIHJrNG0oZG91YmxlICp5LCBkb3VibGUgKmYsIGludCBOLCBkb3VibGUgYSwgZG91YmxlIGIsIGRvdWJsZSBzdGVwLCB2b2lkICgqRikoZG91YmxlLCBkb3VibGUgKiwgZG91YmxlICopKTsKCmludCBtYWluKHZvaWQpIHsKICAgIGludCBOID0gMjsKICAgIGRvdWJsZSAqeSwgKmYsIGEgPSAwLjAsIGIgPSAxLjAsIHN0ZXAgPSAwLjAxOwoKICAgIHkgPSBkdmVjdG9yKDAsIE4gLSAxKTsKICAgIGYgPSBkdmVjdG9yKDAsIE4gLSAxKTsKCiAgICB5WzBdID0gMC4xOyAKICAgIHlbMV0gPSAwLjA7CgogICAgcHJpbnRmKCJ4XHR5XG4iKTsgLy8g44OY44OD44OA44O844KS5Ye65YqbCiAgICByazRtKHksIGYsIE4sIGEsIGIsIHN0ZXAsIEZVTkMpOyAvLyDjg5XjgqHjgqTjg6vmk43kvZzjgpLjgZvjgZrmqJnmupblh7rlipsKCiAgICBmcmVlX2R2ZWN0b3IoeSk7CiAgICBmcmVlX2R2ZWN0b3IoZik7CgogICAgcmV0dXJuIDA7Cn0KCnZvaWQgRlVOQyhkb3VibGUgeCwgZG91YmxlICp5LCBkb3VibGUgKmYpIHsKICAgIGZbMF0gPSB5WzFdOyAvLyBkeTEvZHggPSB5MgogICAgZlsxXSA9IGNvcyh4KSAtIDQgKiB5WzFdIC0gMyAqIHlbMF07IC8vIGR5Mi9keCA9IGNvcyh4KSAtIDR5MiAtIDN5MQp9Cgp2b2lkIHJrNG0oZG91YmxlICp5LCBkb3VibGUgKmYsIGludCBOLCBkb3VibGUgYSwgZG91YmxlIGIsIGRvdWJsZSBzdGVwLCB2b2lkICgqRikoZG91YmxlLCBkb3VibGUgKiwgZG91YmxlICopKSB7CiAgICBkb3VibGUgKmsxLCAqazIsICprMywgKms0LCAqdG1wLCB4LCBoOwogICAgaW50IGo7CgogICAgazEgPSBkdmVjdG9yKDAsIE4gLSAxKTsKICAgIGsyID0gZHZlY3RvcigwLCBOIC0gMSk7CiAgICBrMyA9IGR2ZWN0b3IoMCwgTiAtIDEpOwogICAgazQgPSBkdmVjdG9yKDAsIE4gLSAxKTsKICAgIHRtcCA9IGR2ZWN0b3IoMCwgTiAtIDEpOwoKICAgIGggPSBzdGVwOwogICAgeCA9IGE7CgogICAgd2hpbGUgKHggPD0gYikgewogICAgICAgIHByaW50ZigiJS4ybGZcdCUuMTBsZlxuIiwgeCwgeVswXSk7IC8vIHgg44GoIHlbMF0g44KS5qiZ5rqW5Ye65YqbCgogICAgICAgIEYoeCwgeSwgZik7CiAgICAgICAgZm9yIChqID0gMDsgaiA8IE47IGorKykgazFbal0gPSBmW2pdOwoKICAgICAgICBmb3IgKGogPSAwOyBqIDwgTjsgaisrKSB0bXBbal0gPSB5W2pdICsgaCAqIGsxW2pdIC8gMi4wOwogICAgICAgIEYoeCArIGggLyAyLjAsIHRtcCwgZik7CiAgICAgICAgZm9yIChqID0gMDsgaiA8IE47IGorKykgazJbal0gPSBmW2pdOwoKICAgICAgICBmb3IgKGogPSAwOyBqIDwgTjsgaisrKSB0bXBbal0gPSB5W2pdICsgaCAqIGsyW2pdIC8gMi4wOwogICAgICAgIEYoeCArIGggLyAyLjAsIHRtcCwgZik7CiAgICAgICAgZm9yIChqID0gMDsgaiA8IE47IGorKykgazNbal0gPSBmW2pdOwoKICAgICAgICBmb3IgKGogPSAwOyBqIDwgTjsgaisrKSB0bXBbal0gPSB5W2pdICsgaCAqIGszW2pdOwogICAgICAgIEYoeCArIGgsIHRtcCwgZik7CiAgICAgICAgZm9yIChqID0gMDsgaiA8IE47IGorKykgazRbal0gPSBmW2pdOwoKICAgICAgICBmb3IgKGogPSAwOyBqIDwgTjsgaisrKQogICAgICAgICAgICB5W2pdID0geVtqXSArIGggKiAoazFbal0gKyAyICogazJbal0gKyAyICogazNbal0gKyBrNFtqXSkgLyA2LjA7CgogICAgICAgIHggKz0gaDsgLy8geCDjgpLmm7TmlrAKICAgIH0KCiAgICBmcmVlX2R2ZWN0b3IoazEpOwogICAgZnJlZV9kdmVjdG9yKGsyKTsKICAgIGZyZWVfZHZlY3RvcihrMyk7CiAgICBmcmVlX2R2ZWN0b3IoazQpOwogICAgZnJlZV9kdmVjdG9yKHRtcCk7Cn0KCmRvdWJsZSAqZHZlY3Rvcihsb25nIGksIGxvbmcgaikgewogICAgZG91YmxlICphOwogICAgaWYgKChhID0gbWFsbG9jKChqIC0gaSArIDEpICogc2l6ZW9mKGRvdWJsZSkpKSA9PSBOVUxMKSB7CiAgICAgICAgcHJpbnRmKCLjg6Hjg6Ljg6rjgYznorrkv53jgafjgY3jgb7jgZvjgpNcbiIpOwogICAgICAgIGV4aXQoMSk7CiAgICB9CiAgICByZXR1cm4gYTsKfQoKdm9pZCBmcmVlX2R2ZWN0b3IoZG91YmxlICphKSB7CiAgICBmcmVlKGEpOwp9Cg==