#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; // 終了値を 1.0、刻み幅を 0.01 に設定

    y = dvector(0, N - 1);
    f = dvector(0, N - 1);

    y[0] = 0.1; 
    y[1] = 0.0;

    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\t %.10lf\n", x, y[0], y[1]);

        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);
}