diff --git a/CMakeLists.txt b/CMakeLists.txt
index d9fe3c60a83..2fa06236bb5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -140,7 +140,7 @@ else()
     sundials
     DOWNLOAD_EXTRACT_TIMESTAMP ON
     GIT_REPOSITORY https://github.com/LLNL/sundials
-    GIT_TAG        v6.1.1
+    GIT_TAG        v7.1.1
       # adjust this to the version you need
   )
   FetchContent_GetProperties(sundials)
@@ -150,7 +150,7 @@ else()
   endif()
 endif()
 
-
+set(BOOST_NUMERIC_ODEINT_NO_ADAPTORS ON)
 if (STAN_USE_SYSTEM_BOOST)
   find_package(Boost REQUIRED)
 else()
@@ -162,13 +162,13 @@ else()
   set(BOOST_ENABLE_PYTHON OFF)
   set(BOOST_BUILD_TESTS OFF)
   set(CMAKE_BUILD_TYPE Release)
-  set(BOOST_ENABLE_MPI ON)
+  set(BOOST_ENABLE_MPI ${STAN_MPI})
 
   set(BOOST_INCLUDE_LIBRARIES math numeric/odeint lexical_cast optional random mpi serialization)
   FetchContent_Declare(
     Boost
     DOWNLOAD_EXTRACT_TIMESTAMP ON
-    URL https://github.com/boostorg/boost/releases/download/boost-1.85.0/boost-1.85.0-cmake.tar.xz
+    URL https://github.com/boostorg/boost/releases/download/boost-1.86.0/boost-1.86.0-cmake.tar.xz
   )
   FetchContent_MakeAvailable(Boost)
 endif()
diff --git a/stan/math/rev/functor/algebra_solver_fp.hpp b/stan/math/rev/functor/algebra_solver_fp.hpp
index 871f5f1b09e..850c2c67904 100644
--- a/stan/math/rev/functor/algebra_solver_fp.hpp
+++ b/stan/math/rev/functor/algebra_solver_fp.hpp
@@ -11,7 +11,7 @@
 #include <stan/math/prim/fun/to_array_1d.hpp>
 #include <stan/math/prim/fun/to_vector.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
-#include <sundials/sundials_context.h>
+#include <sundials/sundials_context.hpp>
 #include <kinsol/kinsol.h>
 #include <sunmatrix/sunmatrix_dense.h>
 #include <sunlinsol/sunlinsol_dense.h>
diff --git a/stan/math/rev/functor/cvodes_integrator.hpp b/stan/math/rev/functor/cvodes_integrator.hpp
index b36bae594e0..ae99eb2882e 100644
--- a/stan/math/rev/functor/cvodes_integrator.hpp
+++ b/stan/math/rev/functor/cvodes_integrator.hpp
@@ -8,7 +8,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
-#include <sundials/sundials_context.h>
+#include <sundials/sundials_context.hpp>
 #include <cvodes/cvodes.h>
 #include <nvector/nvector_serial.h>
 #include <sunlinsol/sunlinsol_dense.h>
@@ -66,7 +66,7 @@ class cvodes_integrator {
    * Implements the function of type CVRhsFn which is the user-defined
    * ODE RHS passed to CVODES.
    */
-  static int cv_rhs(realtype t, N_Vector y, N_Vector ydot, void* user_data) {
+  static int cv_rhs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) {
     cvodes_integrator* integrator = static_cast<cvodes_integrator*>(user_data);
     integrator->rhs(t, NV_DATA_S(y), NV_DATA_S(ydot));
     return 0;
@@ -76,7 +76,7 @@ class cvodes_integrator {
    * Implements the function of type CVSensRhsFn which is the
    * RHS of the sensitivity ODE system.
    */
-  static int cv_rhs_sens(int Ns, realtype t, N_Vector y, N_Vector ydot,
+  static int cv_rhs_sens(int Ns, sunrealtype t, N_Vector y, N_Vector ydot,
                          N_Vector* yS, N_Vector* ySdot, void* user_data,
                          N_Vector tmp1, N_Vector tmp2) {
     cvodes_integrator* integrator = static_cast<cvodes_integrator*>(user_data);
@@ -90,7 +90,7 @@ class cvodes_integrator {
    * ode_rhs wrt to the states y. The jacobian is stored in column
    * major format.
    */
-  static int cv_jacobian_states(realtype t, N_Vector y, N_Vector fy,
+  static int cv_jacobian_states(sunrealtype t, N_Vector y, N_Vector fy,
                                 SUNMatrix J, void* user_data, N_Vector tmp1,
                                 N_Vector tmp2, N_Vector tmp3) {
     cvodes_integrator* integrator = static_cast<cvodes_integrator*>(user_data);
diff --git a/stan/math/rev/functor/cvodes_integrator_adjoint.hpp b/stan/math/rev/functor/cvodes_integrator_adjoint.hpp
index a1769f4bdee..adf86937242 100644
--- a/stan/math/rev/functor/cvodes_integrator_adjoint.hpp
+++ b/stan/math/rev/functor/cvodes_integrator_adjoint.hpp
@@ -9,7 +9,7 @@
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/functor/for_each.hpp>
-#include <sundials/sundials_context.h>
+#include <sundials/sundials_context.hpp>
 #include <cvodes/cvodes.h>
 #include <nvector/nvector_serial.h>
 #include <sunmatrix/sunmatrix_dense.h>
@@ -576,7 +576,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
    * Implements the function of type CVRhsFn which is the user-defined
    * ODE RHS passed to CVODES.
    */
-  constexpr static int cv_rhs(realtype t, N_Vector y, N_Vector ydot,
+  constexpr static int cv_rhs(sunrealtype t, N_Vector y, N_Vector ydot,
                               void* user_data) {
     return cast_to_self(user_data)->rhs(t, NV_DATA_S(y), NV_DATA_S(ydot));
   }
@@ -611,7 +611,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
    * Implements the function of type CVRhsFnB which is the
    * RHS of the backward ODE system.
    */
-  constexpr static int cv_rhs_adj(realtype t, N_Vector y, N_Vector yB,
+  constexpr static int cv_rhs_adj(sunrealtype t, N_Vector y, N_Vector yB,
                                   N_Vector yBdot, void* user_data) {
     return cast_to_self(user_data)->rhs_adj(t, y, yB, yBdot);
   }
@@ -654,7 +654,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
    * Implements the function of type CVQuadRhsFnB which is the
    * RHS of the backward ODE system's quadrature.
    */
-  constexpr static int cv_quad_rhs_adj(realtype t, N_Vector y, N_Vector yB,
+  constexpr static int cv_quad_rhs_adj(sunrealtype t, N_Vector y, N_Vector yB,
                                        N_Vector qBdot, void* user_data) {
     return cast_to_self(user_data)->quad_rhs_adj(t, y, yB, qBdot);
   }
@@ -693,7 +693,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
    * ode_rhs wrt to the states y. The jacobian is stored in column
    * major format.
    */
-  constexpr static int cv_jacobian_rhs_states(realtype t, N_Vector y,
+  constexpr static int cv_jacobian_rhs_states(sunrealtype t, N_Vector y,
                                               N_Vector fy, SUNMatrix J,
                                               void* user_data, N_Vector tmp1,
                                               N_Vector tmp2, N_Vector tmp3) {
@@ -722,7 +722,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
    * Implements the CVLsJacFnB function for evaluating the jacobian of
    * the adjoint problem wrt to the backward states.
    */
-  constexpr static int cv_jacobian_rhs_adj_states(realtype t, N_Vector y,
+  constexpr static int cv_jacobian_rhs_adj_states(sunrealtype t, N_Vector y,
                                                   N_Vector yB, N_Vector fyB,
                                                   SUNMatrix J, void* user_data,
                                                   N_Vector tmp1, N_Vector tmp2,
diff --git a/stan/math/rev/functor/idas_service.hpp b/stan/math/rev/functor/idas_service.hpp
index 2055c2e2b91..f12e9146de9 100644
--- a/stan/math/rev/functor/idas_service.hpp
+++ b/stan/math/rev/functor/idas_service.hpp
@@ -9,7 +9,7 @@
 #include <nvector/nvector_serial.h>
 #include <sunmatrix/sunmatrix_dense.h>
 #include <sunlinsol/sunlinsol_dense.h>
-#include <sundials/sundials_context.h>
+#include <sundials/sundials_context.hpp>
 #include <ostream>
 #include <vector>
 #include <algorithm>
@@ -106,8 +106,8 @@ struct idas_service {
           "Failed to allocate yps N_Vectors for sensitivities.");
     }
     for (size_t is = 0; is < ns; ++is) {
-      N_VConst(RCONST(0.0), yys[is]);
-      N_VConst(RCONST(0.0), yps[is]);
+      N_VConst(SUN_RCONST(0.0), yys[is]);
+      N_VConst(SUN_RCONST(0.0), yps[is]);
     }
     set_init_sens(yys, yps, n);
     // std::cout << "Calling IDASensInit with ns=" << ns << " inner ns: " << ns_
diff --git a/stan/math/rev/functor/kinsol_data.hpp b/stan/math/rev/functor/kinsol_data.hpp
index a0847cf9fbd..dbd83152eaf 100644
--- a/stan/math/rev/functor/kinsol_data.hpp
+++ b/stan/math/rev/functor/kinsol_data.hpp
@@ -6,7 +6,7 @@
 #include <stan/math/prim/fun/to_array_1d.hpp>
 #include <stan/math/prim/fun/to_vector.hpp>
 #include <stan/math/prim/functor/apply.hpp>
-#include <sundials/sundials_context.h>
+#include <sundials/sundials_context.hpp>
 #include <kinsol/kinsol.h>
 #include <sunmatrix/sunmatrix_dense.h>
 #include <sunlinsol/sunlinsol_dense.h>
diff --git a/stan/math/rev/functor/kinsol_solve.hpp b/stan/math/rev/functor/kinsol_solve.hpp
index 08e21e73822..24688288dec 100644
--- a/stan/math/rev/functor/kinsol_solve.hpp
+++ b/stan/math/rev/functor/kinsol_solve.hpp
@@ -6,7 +6,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/to_array_1d.hpp>
 #include <stan/math/prim/fun/to_vector.hpp>
-#include <sundials/sundials_context.h>
+#include <sundials/sundials_context.hpp>
 #include <kinsol/kinsol.h>
 #include <sunmatrix/sunmatrix_dense.h>
 #include <sunlinsol/sunlinsol_dense.h>