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

  h = (int(2)-int(1))/N; 
  y = zeros (length(y0), N+1);
  y(:,1) = y0;
  for j = 2:N+1
    tn = int(1) + (j-1/2)*h;
    yn = y(:,j-1);
    yi = yn;
    for k = 1:20
      yn = y(:,j-1) + h * feval (f, tn, 0.5*(y(:,j-1)+yi));
      if abs(yn-yi) < 10^-11
	break;
      end
      yi = yn;
    end
    y(:,j) = yn;
  end

end