We present an adaptive method of approximate integration of convex (as well as concave) functions based on a certain refinement of the celebrated Hermite–Hadamard inequality. Numerical experiments are performed and the role of harmonic numbers is shown.