function ys = ode_rk4(fun, ts, y0, params)

ys = zeros(length(y0), length(ts));
ys(:,1) = y0;
for i = 2:length(ts)
  t = ts(i-1);
  y = ys(:,i-1);
  h = ts(i) - t;
  k1 = fun(t      , y       , params);
  k2 = fun(t + h/2, y+h/2*k1, params);
  k3 = fun(t + h/2, y+h/2*k2, params);
  k4 = fun(t + h  , y+h*k3  , params);
  ys(:,i) = y + h*( k1/6 + k2/3 + k3/3 + k4/6 );
end