function ys = ode_rk2(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;
  ys(:,i) = y + h*fun(t + h/2, y + h/2*fun(t, y, params), params);
end