diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2550a0a05e..341bb29dd2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -269,16 +269,10 @@ target_link_libraries(
          SUNDIALS::sunmatrixband_static
          SUNDIALS::sunmatrixdense_static
          SUNDIALS::sunmatrixsparse_static
-         SUNDIALS::sunlinsolband_static
-         SUNDIALS::sunlinsoldense_static
-         SUNDIALS::sunlinsolpcg_static
          SUNDIALS::sunlinsolspbcgs_static
          SUNDIALS::sunlinsolspfgmr_static
          SUNDIALS::sunlinsolspgmr_static
          SUNDIALS::sunlinsolsptfqmr_static
-         SUNDIALS::sunlinsolklu_static
-         SUNDIALS::sunnonlinsolnewton_static
-         SUNDIALS::sunnonlinsolfixedpoint_static
          SUNDIALS::cvodes_static
          SUNDIALS::idas_static
          $<$<BOOL:${Boost_CHRONO_FOUND}>:Boost::chrono>
@@ -286,6 +280,12 @@ target_link_libraries(
          ${CMAKE_DL_LIBS}
   PRIVATE
     BLAS::BLAS
+    SUNDIALS::sunnonlinsolnewton_static
+    SUNDIALS::sunnonlinsolfixedpoint_static
+    SUNDIALS::sunlinsolklu_static
+    SUNDIALS::sunlinsolband_static
+    SUNDIALS::sunlinsoldense_static
+    SUNDIALS::sunlinsolpcg_static
     $<$<BOOL:${SUNDIALS_SUPERLUMT_ENABLE}>:SUNDIALS::sundials_sunlinsolsuperlumt>
 )
 
diff --git a/include/amici/abstract_model.h b/include/amici/abstract_model.h
index f340214e7e..04208de1d9 100644
--- a/include/amici/abstract_model.h
+++ b/include/amici/abstract_model.h
@@ -6,10 +6,6 @@
 #include "amici/sundials_matrix_wrapper.h"
 #include "amici/vector.h"
 
-#include <sunmatrix/sunmatrix_band.h>
-#include <sunmatrix/sunmatrix_dense.h>
-#include <sunmatrix/sunmatrix_sparse.h>
-
 #include <memory>
 
 namespace amici {
diff --git a/include/amici/forwardproblem.h b/include/amici/forwardproblem.h
index 67658ebee4..ffcca63666 100644
--- a/include/amici/forwardproblem.h
+++ b/include/amici/forwardproblem.h
@@ -5,9 +5,7 @@
 #include "amici/misc.h"
 #include "amici/model.h"
 #include "amici/vector.h"
-#include <amici/amici.h>
 
-#include <sundials/sundials_direct.h>
 #include <vector>
 
 namespace amici {
diff --git a/include/amici/misc.h b/include/amici/misc.h
index 6874c19185..b9f94ece1e 100644
--- a/include/amici/misc.h
+++ b/include/amici/misc.h
@@ -4,7 +4,6 @@
 #include "amici/defines.h"
 #include "amici/exception.h"
 #include "amici/vector.h"
-#include <sunmatrix/sunmatrix_sparse.h> // SUNMatrixContent_Sparse
 
 #include <algorithm>
 #include <ctime>
diff --git a/include/amici/model_dae.h b/include/amici/model_dae.h
index 2e2e7f2bcc..15f96c890e 100644
--- a/include/amici/model_dae.h
+++ b/include/amici/model_dae.h
@@ -3,12 +3,8 @@
 
 #include "amici/model.h"
 
-#include <nvector/nvector_serial.h>
-
-#include <sunmatrix/sunmatrix_band.h>
-#include <sunmatrix/sunmatrix_dense.h>
-#include <sunmatrix/sunmatrix_sparse.h>
-
+#include <sundials/sundials_nvector.h>
+#include <sundials/sundials_matrix.h>
 #include <numeric>
 #include <vector>
 
diff --git a/include/amici/model_ode.h b/include/amici/model_ode.h
index 8386eb929e..4cb71814ef 100644
--- a/include/amici/model_ode.h
+++ b/include/amici/model_ode.h
@@ -3,12 +3,8 @@
 
 #include "amici/model.h"
 
-#include <nvector/nvector_serial.h>
-
+#include <sundials/sundials_nvector.h>
 #include <sundials/sundials_matrix.h>
-#include <sunmatrix/sunmatrix_band.h>
-#include <sunmatrix/sunmatrix_dense.h>
-#include <sunmatrix/sunmatrix_sparse.h>
 
 #include <vector>
 
diff --git a/include/amici/newton_solver.h b/include/amici/newton_solver.h
index 3559abb711..9c450ac97b 100644
--- a/include/amici/newton_solver.h
+++ b/include/amici/newton_solver.h
@@ -1,8 +1,9 @@
 #ifndef amici_newton_solver_h
 #define amici_newton_solver_h
 
-#include "amici/solver.h"
+#include "amici/defines.h"
 #include "amici/vector.h"
+#include "amici/sundials_linsol_wrapper.h"
 
 namespace amici {
 
diff --git a/include/amici/solver_cvodes.h b/include/amici/solver_cvodes.h
index fc91284240..1d6339b0b5 100644
--- a/include/amici/solver_cvodes.h
+++ b/include/amici/solver_cvodes.h
@@ -5,8 +5,6 @@
 #include "amici/solver.h"
 #include "amici/vector.h"
 
-#include <sundials/sundials_matrix.h>
-
 namespace amici {
 class ExpData;
 class ReturnData;
diff --git a/include/amici/solver_idas.h b/include/amici/solver_idas.h
index c10e191f4b..8aad805472 100644
--- a/include/amici/solver_idas.h
+++ b/include/amici/solver_idas.h
@@ -3,8 +3,6 @@
 
 #include "amici/solver.h"
 
-#include <sundials/sundials_matrix.h>
-
 namespace amici {
 class ExpData;
 class ReturnData;
diff --git a/include/amici/steadystateproblem.h b/include/amici/steadystateproblem.h
index f680b6c128..f2b4ccf2bd 100644
--- a/include/amici/steadystateproblem.h
+++ b/include/amici/steadystateproblem.h
@@ -6,8 +6,6 @@
 #include <amici/newton_solver.h>
 #include <amici/vector.h>
 
-#include <nvector/nvector_serial.h>
-
 #include <memory>
 
 namespace amici {
@@ -15,6 +13,7 @@ namespace amici {
 class ExpData;
 class Solver;
 class Model;
+class BackwardProblem;
 
 /**
  * @brief The SteadystateProblem class solves a steady-state problem using
diff --git a/include/amici/sundials_linsol_wrapper.h b/include/amici/sundials_linsol_wrapper.h
index d17633cd6b..dc1d5c7d83 100644
--- a/include/amici/sundials_linsol_wrapper.h
+++ b/include/amici/sundials_linsol_wrapper.h
@@ -5,19 +5,12 @@
 #include "amici/vector.h"
 
 #include <sundials/sundials_config.h>
-#include <sunlinsol/sunlinsol_band.h>
-#include <sunlinsol/sunlinsol_dense.h>
-#include <sunlinsol/sunlinsol_klu.h>
-#include <sunlinsol/sunlinsol_pcg.h>
+#include <sundials/sundials_iterative.h>
 #include <sunlinsol/sunlinsol_spbcgs.h>
 #include <sunlinsol/sunlinsol_spfgmr.h>
 #include <sunlinsol/sunlinsol_spgmr.h>
 #include <sunlinsol/sunlinsol_sptfqmr.h>
-#ifdef SUNDIALS_SUPERLUMT
-#include <sunlinsol/sunlinsol_superlumt.h>
-#endif
-#include <sunnonlinsol/sunnonlinsol_fixedpoint.h>
-#include <sunnonlinsol/sunnonlinsol_newton.h>
+#include <sundials/sundials_nonlinearsolver.h>
 
 namespace amici {
 
diff --git a/src/sundials_linsol_wrapper.cpp b/src/sundials_linsol_wrapper.cpp
index 0f411bad63..f1133855af 100644
--- a/src/sundials_linsol_wrapper.cpp
+++ b/src/sundials_linsol_wrapper.cpp
@@ -1,6 +1,15 @@
 #include <amici/sundials_linsol_wrapper.h>
 
 #include <amici/exception.h>
+#include <sunlinsol/sunlinsol_klu.h>
+#include <sunlinsol/sunlinsol_band.h>
+#include <sunlinsol/sunlinsol_dense.h>
+#include <sunlinsol/sunlinsol_pcg.h>
+#ifdef SUNDIALS_SUPERLUMT
+#include <sunlinsol/sunlinsol_superlumt.h>
+#endif
+#include <sunnonlinsol/sunnonlinsol_fixedpoint.h>
+#include <sunnonlinsol/sunnonlinsol_newton.h>
 
 #include <utility>