solve_ivp defaults are rtol=1e-3: 28% error after 160 periods; failures don't raise
SciPy 1.18.1, harmonic oscillator y''=-y over t in [0, 1000]. The default RK45 (rtol=1e-3, atol=1e-6) ends with |y(T)-cos(T)| = 0.278. rtol=1e-6 gives 4.0e-4. rtol=1e-10, atol=1e-12 gives 7.4e-9 with 191k evaluations; DOP853 at the same tolerances gives 5.6e-9 with 37k. odeint with only [0, T] in t gives error 0.31 and an 'Excess work done' warning, because it is limited by mxstep, not accuracy. On a stiff problem (Robertson, t to 100) RK45 stops at t=1.37 with status -1 ('Required step size is less than spacing between numbers'). It does not raise, so always check sol.success; BDF/Radau/LSODA finish with 123/208/134 evaluations. Without t_eval or dense_output the default output is only the accepted steps (15 points on [0, 10]), so plots look jagged. That is not an accuracy problem.
Every number reproduces with SciPy 1.18.1 / NumPy 2.5.4 on Python 3.12.16: errors 0.28, 4.0e-4, 7.4e-9 (191234 evaluations) and 5.6e-9 (DOP853, 37418), odeint 0.31 with 'Excess work done', RK45 on Robertson stopping at t = 1.37 with status -1, BDF/Radau/LSODA at 123/208/134 evaluations, and 15 output points on [0, 10]. I also checked the mxstep explanation. With mxstep=5000 odeint still warns. With mxstep=100000 it finishes cleanly with error 1.3e-5 (14051 evaluations). Passing a dense t (np.linspace(0, 1000, 2001)) at the default mxstep gives the same 1.3e-5 with no warning, since mxstep is a limit per output point in t (per the docstring). So a finer t grid fixes it as well as a larger mxstep does.