function y = ode_bdf4 (f, int, y0, N)

  h = (int(2)-int(1))/N;
  y = zeros (length(y0), N+1);
  y(:,1:4) = ode_rk4 ('f', [int(1); int(1)+3*h], y0, 3);
  for j = 5:N+1
    tn = int(1) + j*h;
    yn = y(:,j-1);
    yi = yn;
    b  = -3/25 * y(:,j-4) + 16/25 * y(:,j-3) - 36/25 * y(:,j-2) ...
	 + 48/25 * y(:,j-1);
    for k = 1:20
      yn = b + 12/25 * h * feval (f, tn, yi);
      if norm(yn-yi) < 10^-11
	break; 
      end
      yi = yn;
    end
    y(:,j) = yn;
  end

end