function y = ode_bdf4_broyden (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);
  B = eye(length(y0)); % approximate Jacobian
  for j = 5:N+1
    tn = int(1) + j*h;
    b  = -3/25 * y(:,j-4) + 16/25 * y(:,j-3) - 36/25 * y(:,j-2) ...
	 + 48/25 * y(:,j-1);
    yi = y(:,j-1);
    fi = b + 12/25 * h * feval (f, tn, yi) - yi;
    for k = 1:20
      dy = -B\fi;
      yi = yi + dy;
      if norm(dy) < 10^-11
	break;
      end
      fj = b + 12/25 * h * feval (f, tn, yi) - yi;
      df = fj - fi;
      fi = fj;
      B = B + ((df - B*dy) * dy')/(dy' * dy);     
    end
    y(:,j) = yi;
  end

end