#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
// ベクトル領域を確保する関数
// i: 開始インデックス, j: 終了インデックス
double *dvector(long i, long j) {
    double *a;
    if ((a = malloc(((j - i + 1) * sizeof(double)))) == NULL) { // メモリ確保
        printf("メモリが確保できません\n");
        exit(1);
    }
    return (a - i); // インデックスを 1 ベースで使えるように調整
}
 
// 確保したベクトル領域を解放する関数
void free_dvector(double *a, long i) {
    free((void *)(a + i));
}
 
// 連立微分方程式の右辺を定義する関数
// x: 現在の x 値, y: 現在の y ベクトル, f: 微分方程式の結果を格納
void FUNC(double x, double *y, double *f) {
    f[1] = y[2] - 3 * y[1]; 
    f[2] = -y[2] + cos(x);  
}
 
// 4次のルンゲ＝クッタ法
// y: 解ベクトル, f: 傾きベクトル, N: ベクトルの次元, a: 開始点, b: 終了点, n: 分割数, F: 右辺関数
void rk4m(double *y, double *f, int N, double a, double b, int n, void (*F)()) {
    double *k1, *k2, *k3, *k4, *tmp, x, h;
    int i, j;
 
    // ベクトル領域の確保
    k1 = dvector(1, N);
    k2 = dvector(1, N);
    k3 = dvector(1, N);
    k4 = dvector(1, N);
    tmp = dvector(1, N);
 
    h = (b - a) / n; // ステップ幅
    x = a; // x の初期値
    printf("%.3lf\t %.3lf\t %.10lf\n", x, y[1], y[2]); // 初期条件を出力
 
    // ループで数値解を計算
    for (i = 0; i < n; i++) {
        // k1 の計算
        (*F)(x, y, f); // f = F(x, y)
        for (j = 1; j <= N; j++) k1[j] = f[j];
 
        // k2 の計算
        for (j = 1; j <= N; j++) tmp[j] = y[j] + h * k1[j] / 2.0;
        (*F)(x + h / 2.0, tmp, f);
        for (j = 1; j <= N; j++) k2[j] = f[j];
 
        // k3 の計算
        for (j = 1; j <= N; j++) tmp[j] = y[j] + h * k2[j] / 2.0;
        (*F)(x + h / 2.0, tmp, f);
        for (j = 1; j <= N; j++) k3[j] = f[j];
 
        // k4 の計算
        for (j = 1; j <= N; j++) tmp[j] = y[j] + h * k3[j];
        (*F)(x + h, tmp, f);
        for (j = 1; j <= N; j++) k4[j] = f[j];
 
        // 次の y の計算
        for (j = 1; j <= N; j++) {
            y[j] = y[j] + h * (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6.0;
        }
 
        x += h; // 次の x 値
        printf("%.3lf\t %.10lf\t %.10lf\n", x, y[1], y[2]); // 結果を出力
    }
 
    // ベクトル領域の解放
    free_dvector(k1, 1);
    free_dvector(k2, 1);
    free_dvector(k3, 1);
    free_dvector(k4, 1);
    free_dvector(tmp, 1);
}
 
// メイン関数
int main(void) {
    int n, N = 2; // n: 分割数, N: 連立方程式の次元
    double *y, *f, a = 0.0, b = 5.0; // a: 開始点, b: 終了点
 
    y = dvector(1, N); // 解ベクトルの確保
    f = dvector(1, N); // 傾きベクトルの確保
 
    // 初期条件
    y[1] = 0.1; // y1(初期値)
    y[2] = 0.0; // y2(初期値)
 
    printf("分割数を入力してください\n");
    scanf("%d", &n); // 分割数の入力
 
    // ルンゲ＝クッタ法を実行
    rk4m(y, f, N, a, b, n, FUNC);
 
    // メモリ解放
    free_dvector(y, 1);
    free_dvector(f, 1);
 
    return 0;
}
 

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KI2luY2x1ZGUgPG1hdGguaD4KCi8vIOODmeOCr+ODiOODq+mgmOWfn+OCkueiuuS/neOBmeOCi+mWouaVsAovLyBpOiDplovlp4vjgqTjg7Pjg4fjg4Pjgq/jgrksIGo6IOe1guS6huOCpOODs+ODh+ODg+OCr+OCuQpkb3VibGUgKmR2ZWN0b3IobG9uZyBpLCBsb25nIGopIHsKICAgIGRvdWJsZSAqYTsKICAgIGlmICgoYSA9IG1hbGxvYygoKGogLSBpICsgMSkgKiBzaXplb2YoZG91YmxlKSkpKSA9PSBOVUxMKSB7IC8vIOODoeODouODqueiuuS/nQogICAgICAgIHByaW50Zigi44Oh44Oi44Oq44GM56K65L+d44Gn44GN44G+44Gb44KTXG4iKTsKICAgICAgICBleGl0KDEpOwogICAgfQogICAgcmV0dXJuIChhIC0gaSk7IC8vIOOCpOODs+ODh+ODg+OCr+OCueOCkiAxIOODmeODvOOCueOBp+S9v+OBiOOCi+OCiOOBhuOBq+iqv+aVtAp9CgovLyDnorrkv53jgZfjgZ/jg5njgq/jg4jjg6vpoJjln5/jgpLop6PmlL7jgZnjgovplqLmlbAKdm9pZCBmcmVlX2R2ZWN0b3IoZG91YmxlICphLCBsb25nIGkpIHsKICAgIGZyZWUoKHZvaWQgKikoYSArIGkpKTsKfQoKLy8g6YCj56uL5b6u5YiG5pa556iL5byP44Gu5Y+z6L6644KS5a6a576p44GZ44KL6Zai5pWwCi8vIHg6IOePvuWcqOOBriB4IOWApCwgeTog54++5Zyo44GuIHkg44OZ44Kv44OI44OrLCBmOiDlvq7liIbmlrnnqIvlvI/jga7ntZDmnpzjgpLmoLzntI0Kdm9pZCBGVU5DKGRvdWJsZSB4LCBkb3VibGUgKnksIGRvdWJsZSAqZikgewogICAgZlsxXSA9IHlbMl0gLSAzICogeVsxXTsgCiAgICBmWzJdID0gLXlbMl0gKyBjb3MoeCk7ICAKfQoKLy8gNOasoeOBruODq+ODs+OCsu+8neOCr+ODg+OCv+azlQovLyB5OiDop6Pjg5njgq/jg4jjg6ssIGY6IOWCvuOBjeODmeOCr+ODiOODqywgTjog44OZ44Kv44OI44Or44Gu5qyh5YWDLCBhOiDplovlp4vngrksIGI6IOe1guS6hueCuSwgbjog5YiG5Ymy5pWwLCBGOiDlj7PovrrplqLmlbAKdm9pZCByazRtKGRvdWJsZSAqeSwgZG91YmxlICpmLCBpbnQgTiwgZG91YmxlIGEsIGRvdWJsZSBiLCBpbnQgbiwgdm9pZCAoKkYpKCkpIHsKICAgIGRvdWJsZSAqazEsICprMiwgKmszLCAqazQsICp0bXAsIHgsIGg7CiAgICBpbnQgaSwgajsKCiAgICAvLyDjg5njgq/jg4jjg6vpoJjln5/jga7norrkv50KICAgIGsxID0gZHZlY3RvcigxLCBOKTsKICAgIGsyID0gZHZlY3RvcigxLCBOKTsKICAgIGszID0gZHZlY3RvcigxLCBOKTsKICAgIGs0ID0gZHZlY3RvcigxLCBOKTsKICAgIHRtcCA9IGR2ZWN0b3IoMSwgTik7CgogICAgaCA9IChiIC0gYSkgLyBuOyAvLyDjgrnjg4bjg4Pjg5fluYUKICAgIHggPSBhOyAvLyB4IOOBruWIneacn+WApAogICAgcHJpbnRmKCIlLjNsZlx0ICUuM2xmXHQgJS4xMGxmXG4iLCB4LCB5WzFdLCB5WzJdKTsgLy8g5Yid5pyf5p2h5Lu244KS5Ye65YqbCgogICAgLy8g44Or44O844OX44Gn5pWw5YCk6Kej44KS6KiI566XCiAgICBmb3IgKGkgPSAwOyBpIDwgbjsgaSsrKSB7CiAgICAgICAgLy8gazEg44Gu6KiI566XCiAgICAgICAgKCpGKSh4LCB5LCBmKTsgLy8gZiA9IEYoeCwgeSkKICAgICAgICBmb3IgKGogPSAxOyBqIDw9IE47IGorKykgazFbal0gPSBmW2pdOwoKICAgICAgICAvLyBrMiDjga7oqIjnrpcKICAgICAgICBmb3IgKGogPSAxOyBqIDw9IE47IGorKykgdG1wW2pdID0geVtqXSArIGggKiBrMVtqXSAvIDIuMDsKICAgICAgICAoKkYpKHggKyBoIC8gMi4wLCB0bXAsIGYpOwogICAgICAgIGZvciAoaiA9IDE7IGogPD0gTjsgaisrKSBrMltqXSA9IGZbal07CgogICAgICAgIC8vIGszIOOBruioiOeulwogICAgICAgIGZvciAoaiA9IDE7IGogPD0gTjsgaisrKSB0bXBbal0gPSB5W2pdICsgaCAqIGsyW2pdIC8gMi4wOwogICAgICAgICgqRikoeCArIGggLyAyLjAsIHRtcCwgZik7CiAgICAgICAgZm9yIChqID0gMTsgaiA8PSBOOyBqKyspIGszW2pdID0gZltqXTsKCiAgICAgICAgLy8gazQg44Gu6KiI566XCiAgICAgICAgZm9yIChqID0gMTsgaiA8PSBOOyBqKyspIHRtcFtqXSA9IHlbal0gKyBoICogazNbal07CiAgICAgICAgKCpGKSh4ICsgaCwgdG1wLCBmKTsKICAgICAgICBmb3IgKGogPSAxOyBqIDw9IE47IGorKykgazRbal0gPSBmW2pdOwoKICAgICAgICAvLyDmrKHjga4geSDjga7oqIjnrpcKICAgICAgICBmb3IgKGogPSAxOyBqIDw9IE47IGorKykgewogICAgICAgICAgICB5W2pdID0geVtqXSArIGggKiAoazFbal0gKyAyICogazJbal0gKyAyICogazNbal0gKyBrNFtqXSkgLyA2LjA7CiAgICAgICAgfQoKICAgICAgICB4ICs9IGg7IC8vIOasoeOBriB4IOWApAogICAgICAgIHByaW50ZigiJS4zbGZcdCAlLjEwbGZcdCAlLjEwbGZcbiIsIHgsIHlbMV0sIHlbMl0pOyAvLyDntZDmnpzjgpLlh7rlipsKICAgIH0KCiAgICAvLyDjg5njgq/jg4jjg6vpoJjln5/jga7op6PmlL4KICAgIGZyZWVfZHZlY3RvcihrMSwgMSk7CiAgICBmcmVlX2R2ZWN0b3IoazIsIDEpOwogICAgZnJlZV9kdmVjdG9yKGszLCAxKTsKICAgIGZyZWVfZHZlY3RvcihrNCwgMSk7CiAgICBmcmVlX2R2ZWN0b3IodG1wLCAxKTsKfQoKLy8g44Oh44Kk44Oz6Zai5pWwCmludCBtYWluKHZvaWQpIHsKICAgIGludCBuLCBOID0gMjsgLy8gbjog5YiG5Ymy5pWwLCBOOiDpgKPnq4vmlrnnqIvlvI/jga7mrKHlhYMKICAgIGRvdWJsZSAqeSwgKmYsIGEgPSAwLjAsIGIgPSA1LjA7IC8vIGE6IOmWi+Wni+eCuSwgYjog57WC5LqG54K5CgogICAgeSA9IGR2ZWN0b3IoMSwgTik7IC8vIOino+ODmeOCr+ODiOODq+OBrueiuuS/nQogICAgZiA9IGR2ZWN0b3IoMSwgTik7IC8vIOWCvuOBjeODmeOCr+ODiOODq+OBrueiuuS/nQoKICAgIC8vIOWIneacn+adoeS7tgogICAgeVsxXSA9IDAuMTsgLy8geTEo5Yid5pyf5YCkKQogICAgeVsyXSA9IDAuMDsgLy8geTIo5Yid5pyf5YCkKQoKICAgIHByaW50Zigi5YiG5Ymy5pWw44KS5YWl5Yqb44GX44Gm44GP44Gg44GV44GEXG4iKTsKICAgIHNjYW5mKCIlZCIsICZuKTsgLy8g5YiG5Ymy5pWw44Gu5YWl5YqbCgogICAgLy8g44Or44Oz44Ky77yd44Kv44OD44K/5rOV44KS5a6f6KGMCiAgICByazRtKHksIGYsIE4sIGEsIGIsIG4sIEZVTkMpOwoKICAgIC8vIOODoeODouODquino+aUvgogICAgZnJlZV9kdmVjdG9yKHksIDEpOwogICAgZnJlZV9kdmVjdG9yKGYsIDEpOwoKICAgIHJldHVybiAwOwp9Cg==

NTA=

50