% Takes one time step of size dt using the rungekutta method
% starting from v_0 and where the function F(v,t) gives the
% time derivatives.
function v = rungekutta_6(v, t , dt, F)
    s = 7
    k = zeros(length(v),s)
    a = zeros(7,6);
    c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1];
    b = [1/12, 0, 0, 0, 5/12, 5/12, 1/12];
    a = [
        0,                           0,                          0,                       0,                     0,                 0;
        4/7,                         0,                          0,                       0,                     0,                 0;
        115/112,                     -5/16,                      0,                       0,                     0,                 0;
        589/630,                     5/18,                       -16/45,                  0,                     0,                 0;
        229/1200 - 29/6000*sqrt(5),  119/240 - 187/1200*sqrt(5), -14/75 + 34/375*sqrt(5), -3/100*sqrt(5),        0,                 0;
        71/2400 - 587/12000*sqrt(5), 187/480 - 391/2400*sqrt(5), -38/75 + 26/375*sqrt(5), 27/80 - 3/400*sqrt(5), (1+sqrt(5))/4,     0;
        -49/480 + 43/160*sqrt(5),    -425/96 + 51/32*sqrt(5),    52/15 - 4/5*sqrt(5),     -27/16 + 3/16*sqrt(5), 5/4 - 3/4*sqrt(5), 5/2 - 1/2*sqrt(5);
    ];

    for i = 1:s
        u = v;
        for j = 1:i-1
            u = u + dt*a(i,j)*k(:,j);
        end
        k(:,i) = F(t+c(i)*dt,u);
    end

    for i = 1:s
        v = v + dt*b(i)*k(:,i);
    end
end