RK

// Compile: gcc ergasia1_RKF4_2nd_order.c -lm -o run

#include <stdio.h>
#include <math.h> // add -lm flag at compilation 

#define max(a, b) a > b ? a : b
#define min(a, b) a < b ? a : b


double func(double t, double y) {
    return 2 * (y - 1) * (y - 1); // Set the function of your preference
}


int main() {
    double t0, tF, y0, y0_deriv;
    double Rmax = 0.0001, ScaleMin = 0.1, ScaleMax = 4.0; // default values.
    
    printf("Enter t0, tF, y(0), y'(0): ");
    scanf("%lf%lf%lf", &t0, &tF, &y0, &y0_deriv);

    double h = pow(Rmax, 1.0 / 4.0);
    double Hmin = h * 0.0001;
    double t = t0, y = y0;

    do {
        if (t + h > tF) {
            h = tF - t;
        }

        double m1 = h * func(t, y);
        double m2 = h * func(t + 1.0 / 4.0 * h, y + 1.0 / 4.0 * m1);
        double m3 = h * func(t + 3.0 / 8.0 * h, y + 3.0 / 32.0 * m1 + 9.0 / 32.0 * m2);
        double m4 = h * func(t + 12.0 / 13.0 * h, y + 1932.0 / 2197.0 * m1 - 7200.0 / 2197.0 * m2 + 7296.0 / 2197.0 * m3);
        double m5 = h * func(t + h, y + 439.0 / 216.0 * m1 - 8.0 * m2 + 3680.0 / 513.0 * m3 - 845.0 / 4104.0 * m4);
        double m6 = h * func(t + 1.0 / 2.0 * h, y - 8.0 / 27.0 * m1 + 2.0 * m2 - 3544.0 / 2565.0 * m3 + 1859.0 / 4104.0 * m4 - 11.0 / 40.0 * m5);
        double error_estimate = 1.0 / 360.0 * m1 - 128.0 / 4275.0 * m2 - 2197.0 / 75240.0 * m4 + 1.0 / 50.0 * m5 + 2.0 / 55.0 * m6;

        double ratio = fabs(error_estimate) / h;

        if (ratio < Rmax) {
            t = t + h;
            y += 25.0 / 216.0 * m1 + 1408.0 / 2565.0 * m3 + 2197.0 / 4104.0 * m4 - 1.0 / 5.0 * m5;

            printf("%.10lf %.10lf %.10lf\n", t, h, y);
        }
        
        double ScaleFactor = 0.84 * pow(Rmax / ratio, 1.0 / 4.0);

        ScaleFactor = max(ScaleFactor, ScaleMin);
        ScaleFactor = min(ScaleFactor, ScaleMax);

        h *= ScaleFactor;
    } while (t != tF && h >= Hmin);

    if (t == tF) {
        printf("%.10lf %.10lf\n", t, y);
    } else {
        printf("συµβαίνει να είναι h < Hmin, κάποια ιδιαιτερότητα παρατηρείται ”κοντά” στο τρέχον t\n");
    }

    return 0;
}