function [T,X] = rungekutta45(fun, tlims, x0, dt)
N = floor((tlims(2)-tlims(1))/dt) + 1;
T = zeros(N,1);
X = zeros(N,1);
T(1) = tlims(1);
X(1) = x0;
for n = 1:N
    T(n+1) = T(n) + dt;
    K1 = feval(fun,X(n))*dt;
    K2 = feval(fun,X(n)+K1/2)*dt;
    K3 = feval(fun,X(n)+K2/2)*dt;
    K4 = feval(fun,X(n)+K3)*dt;
    X(n+1) = X(n) + (K1 + 2*K2 + 2*K3 + K4)/6;
end