fn runge_kutta4(fx: &Fn(f64, f64) -> f64, x: f64, y: f64, dx: f64) -> f64 {
    let k1 = dx * fx(x, y);
    let k2 = dx * fx(x + dx / 2.0, y + k1 / 2.0);
    let k3 = dx * fx(x + dx / 2.0, y + k2 / 2.0);
    let k4 = dx * fx(x + dx, y + k3);

    y + (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0
}

fn f(x: f64, y: f64) -> f64 {
    x * y.sqrt()
}

fn actual(x: f64) -> f64 {
    (1.0 / 16.0) * (x * x + 4.0).powi(2)
}

fn main() {
    let mut y = 1.0;
    let mut x = 0.0;
    let step = 0.1;
    let max_steps = 101;
    let sample_every_n = 10;

    for steps in 0..max_steps {
        if steps % sample_every_n == 0 {
            println!("y({}):\t{:.10}\t\t {:E}", x, y, actual(x) - y)
        }

        y = runge_kutta4(&f, x, y, step);

        x = ((x * 10.0) + (step * 10.0)) / 10.0;
    }
}