function [T U] = odeRK4(func, t0, t_final, u0, N)
% Dag Lindbo, 2011-02-23
    
    T = zeros(N+1,1);
    U = zeros(N+1,length(u0));
    
    t = t0;
    u = u0(:);
    h = (t_final-t0)/N;
    
    T(1) = t;
    U(1,:) = u;
    
    for i=1:N
        k1 = h*func(t    ,u     );
        k2 = h*func(t+h/2,u+k1/2);
        k3 = h*func(t+h/2,u+k2/2);
        k4 = h*func(t+h  ,u+  k3);
        
        u = u + (k1+2*k2+2*k3+k4)/6;
        t = t+h;
        
        T(i+1) = t;
        U(i+1,:) = u';
    end
    