In the context of multiphysics simulations, partitioned methods enable the reuse of existing solvers while preserving modularity. However, achieving high-order accuracy in time and providing adaptive time stepping remains a challenge, especially when dealing with strong coupling conditions. We address these issues by developing and analyzing a high-order multistep coupling scheme tailored for multiphysics applications. Contrary to classical staggered coupling schemes, where the coupling terms are held constant between two successive coupling time points, this new technique uses high-order polynomials in time to predict the evolution of the coupling conditions during a coupling time step. An explicit and an implicit variant arise naturally, and error estimates can be built to dynamically drive the coupling time step. We first present the numerical analysis of convergence and stability of the multistep coupling scheme. We then build an equivalent of the classical Dahlquist test equation for absolute stability analysis of ODEs, but for coupled systems of ODEs. We suggest a way to connect actual coupled PDEs to this simple model equation. The theoretical approach is then assessed on a conjugate heat transfer (CHT) test-case. We then present a first distributed-memory implementation of this multistep scheme within the C++ coupling library \textsc{Cwipi}, targeting a multi-dimensional solid - fluid CHT case. Considerations on the use of the implicit variant of the scheme are also presented, showing that it can be a viable option for strongly coupled scenarios. Finally, various error estimates are compared in terms of their time step adaptation efficiency. This work provides an interesting computational framework to improve existing code-coupled multiphysics simulations, increasing their accuracy and stability.