function y = ladder(ef, y0, h, a, b)
y(1) = y0;
n = (b - a) / h;
x = a : h : b;
fn0 = y0;
fn1 = y0;
fn2 = y0;
fn3 = y0;
for i = 1 : n
  y(i+1) = y(i) + h / 24 * (55 * fn0 - 59 * fn1 + 37 * fn2 - 9 * fn3);
  y(i+1) = y(i) + h / 2 * (feval(ef, x(i), y(i)) + feval(ef, x(i+1), y(i+1)));
  fn3 = fn2;
  fn2 = fn1;
  fn1 = fn0;
  fn0 = y(i+1);
end