pub fn rk4_nd(f: Vec<impl Fn(f64, Vec<f64>) -> f64>, n: u32, alpha: Vec<f64>, a: f64, b: f64) -> Vec<(f64, Vec<f64>)> {
    let m = f.len();

    // Step 1 in Burden
    let h = (b - a) / (n as f64);
    let mut t = a;

    // Step 2 in Burden
    let mut w = alpha.clone();

    // Step 3 in Burden
    let mut output = vec![];
    output.push((t, w.clone()));

    // Step 4 in Burden
    for _ in 0..n {
        // Step 5 in Burden
        let mut k1 = vec![];
        for j in 0..m {
            k1.push(h * f[j](t, w.clone()));
        }

        // Step 6 in Burden
        let mut k2 = vec![];
        let temp: Vec<f64> = k1.clone().iter().zip(&w).map(|(k, w)| w + k / 2.0).collect();
        for j in 0..m {
            k2.push(h * f[j](t + h / 2.0, temp.clone()));
        }

        // Step 7 in Burden
        let mut k3 = vec![];
        let temp: Vec<f64> = k2.clone().iter().zip(&w).map(|(k, w)| w + k / 2.0).collect();
        for j in 0..m {
            k3.push(h * f[j](t + h / 2.0, temp.clone()));
        }

        // Step 8 in Burden
        let mut k4 = vec![];
        let temp: Vec<f64> = k3.clone().iter().zip(&w).map(|(k, w)| w + k).collect();
        for j in 0..m {
            k4.push(h * f[j](t + h, temp.clone()));
        }

        // Step 9 in Burden
        for j in 0..m {
            w[j] += (k1[j] + 2.0 * k2[j] + 2.0 * k3[j] + k4[j]) / 6.0;
        }

        // Step 10 in Burden
        t += h;

        // Step 11 in Burden
        output.push((t, w.clone()));
    }
    output
}

fn f1(t: f64, w: Vec<f64>) -> f64 { -4.0 * w[0] + 3.0 * w[1] + 6.0 }
fn f2(t: f64, w: Vec<f64>) -> f64 { -2.4 * w[0] + 1.6 * w[1] + 3.6 }

fn main() {
    // Burden, Illustration, pg. 331.
    let mut f = vec![];
    f.push(f1);
    f.push(f2);
    let n = 10;
    let a = 0.0;
    let b = 1.0;
    let alpha = vec![0.0, 0.0];

    let result = rk4_nd(f, n, alpha, a, b);
    println!("{:?}", result);
}