diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4902462052..8165a8187d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -14,6 +14,8 @@ BREAKING CHANGES
 * `ReturnDataView.posteq_status` and `ReturnDataView.preeq_status` now
   return `list[SteadyStateStatus]` instead of an `ndarray[int]` of shape
   `(1, 3)`.
+* `AmiVector::getLength` and `AmiVectorArray::getLength` have been renamed to
+  `size` to be more consistent with STL containers.
 * The following deprecated functionality has been removed:
   * The complete MATLAB interface has been removed.
   * `NonlinearSolverIteration::functional` has been removed,
diff --git a/include/amici/vector.h b/include/amici/vector.h
index 4bc2219555..997084ac29 100644
--- a/include/amici/vector.h
+++ b/include/amici/vector.h
@@ -199,7 +199,7 @@ class AmiVector {
      * @brief returns the length of the vector
      * @return length
      */
-    [[nodiscard]] int getLength() const;
+    [[nodiscard]] int size() const;
 
     /**
      * @brief fills vector with zero values
@@ -310,7 +310,7 @@ class AmiVector {
  */
 inline std::ostream& operator<<(std::ostream& os, AmiVector const& v) {
     os << "[";
-    for (int i = 0; i < v.getLength(); ++i) {
+    for (int i = 0; i < v.size(); ++i) {
         if (i > 0)
             os << ", ";
         os << v.at(i);
@@ -429,7 +429,7 @@ class AmiVectorArray {
      * @brief length of AmiVectorArray
      * @return length
      */
-    int getLength() const;
+    int size() const;
 
     /**
      * @brief set every AmiVector in AmiVectorArray to zero
@@ -482,7 +482,7 @@ class AmiVectorArray {
  */
 inline std::ostream& operator<<(std::ostream& os, AmiVectorArray const& arr) {
     os << "[";
-    for (int i = 0; i < arr.getLength(); ++i) {
+    for (int i = 0; i < arr.size(); ++i) {
         if (i > 0)
             os << ", ";
         os << arr[i];
diff --git a/src/model.cpp b/src/model.cpp
index f154d82c60..82c149515d 100644
--- a/src/model.cpp
+++ b/src/model.cpp
@@ -3038,7 +3038,7 @@ void Model::set_steadystate_mask(std::vector<realtype> const& mask) {
 
 const_N_Vector Model::computeX_pos(const_N_Vector x) {
     if (any_state_non_negative_) {
-        for (int ix = 0; ix < derived_state_.x_pos_tmp_.getLength(); ++ix) {
+        for (int ix = 0; ix < derived_state_.x_pos_tmp_.size(); ++ix) {
             derived_state_.x_pos_tmp_.at(ix)
                 = (state_is_non_negative_.at(ix) && NV_Ith_S(x, ix) < 0)
                       ? 0
diff --git a/src/rdata.cpp b/src/rdata.cpp
index 22e9ef37ac..5d88611e11 100644
--- a/src/rdata.cpp
+++ b/src/rdata.cpp
@@ -290,7 +290,7 @@ void ReturnData::processForwardProblem(
         initializeObjectiveFunction(model.hasQuadraticLLH());
 
     auto const& initialState = fwd.getInitialSimulationState();
-    if (initialState.sol.x.getLength() == 0 && model.nx_solver > 0)
+    if (initialState.sol.x.size() == 0 && model.nx_solver > 0)
         return; // if x wasn't set forward problem failed during initialization
 
     model.setModelState(initialState.mod);
diff --git a/src/solver.cpp b/src/solver.cpp
index 0ccf8dd0a1..32436d02bd 100644
--- a/src/solver.cpp
+++ b/src/solver.cpp
@@ -713,8 +713,8 @@ void Solver::applySensitivityTolerances() const {
 }
 
 void Solver::apply_constraints() const {
-    if (constraints_.getLength() != 0
-        && gsl::narrow<int>(constraints_.getLength()) != nx()) {
+    if (constraints_.size() != 0
+        && gsl::narrow<int>(constraints_.size()) != nx()) {
         throw std::invalid_argument(
             "Constraints must have the same size as the state vector."
         );
@@ -1230,11 +1230,11 @@ void Solver::setErrHandlerFn() const {
         );
 }
 
-int Solver::nplist() const { return sx_.getLength(); }
+int Solver::nplist() const { return sx_.size(); }
 
-int Solver::nx() const { return x_.getLength(); }
+int Solver::nx() const { return x_.size(); }
 
-int Solver::nquad() const { return xQB_.getLength(); }
+int Solver::nquad() const { return xQB_.size(); }
 
 bool Solver::getInitDone() const { return initialized_; }
 
diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp
index 89eb45b48b..baa489baeb 100644
--- a/src/solver_cvodes.cpp
+++ b/src/solver_cvodes.cpp
@@ -299,7 +299,7 @@ void CVodeSolver::apply_constraints() const {
 
     int status = CVodeSetConstraints(
         solver_memory_.get(),
-        constraints_.getLength() > 0 ? constraints_.getNVector() : nullptr
+        constraints_.size() > 0 ? constraints_.getNVector() : nullptr
     );
     if (status != CV_SUCCESS) {
         throw CvodeException(status, "CVodeSetConstraints");
diff --git a/src/solver_idas.cpp b/src/solver_idas.cpp
index 8db9431eeb..616f2df829 100644
--- a/src/solver_idas.cpp
+++ b/src/solver_idas.cpp
@@ -301,7 +301,7 @@ void IDASolver::apply_constraints() const {
 
     int status = IDASetConstraints(
         solver_memory_.get(),
-        constraints_.getLength() > 0 ? constraints_.getNVector() : nullptr
+        constraints_.size() > 0 ? constraints_.getNVector() : nullptr
     );
     if (status != IDA_SUCCESS) {
         throw IDAException(status, "IDASetConstraints");
diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp
index bc63e93006..5b320a9725 100644
--- a/src/steadystateproblem.cpp
+++ b/src/steadystateproblem.cpp
@@ -40,7 +40,7 @@ realtype getWrmsNorm(
     N_VInv(ewt.getNVector(), ewt.getNVector());
 
     // wrms = sqrt(sum((xdot/ewt)**2)/n) where n = size of state vector
-    if (mask.getLength()) {
+    if (mask.size()) {
         return N_VWrmsNormMask(
             const_cast<N_Vector>(xdot.getNVector()), ewt.getNVector(),
             const_cast<N_Vector>(mask.getNVector())
@@ -189,7 +189,7 @@ bool NewtonsMethod::has_converged(
     auto nonnegative = model_->getStateIsNonNegative();
     Expects(nonnegative.size() == state.x.getVector().size());
     auto state_modified = false;
-    for (int ix = 0; ix < state.x.getLength(); ix++) {
+    for (int ix = 0; ix < state.x.size(); ix++) {
         if (state.x[ix] < 0.0 && nonnegative[ix]) {
             state.x[ix] = 0.0;
             state_modified = true;
diff --git a/src/sundials_linsol_wrapper.cpp b/src/sundials_linsol_wrapper.cpp
index 96dd480515..c3eb1fbba9 100644
--- a/src/sundials_linsol_wrapper.cpp
+++ b/src/sundials_linsol_wrapper.cpp
@@ -181,7 +181,7 @@ SUNLinSolBand::SUNLinSolBand(N_Vector x, SUNMatrixWrapper A)
 
 SUNLinSolBand::SUNLinSolBand(AmiVector const& x, int ubw, int lbw)
     : SUNLinSolWrapper(
-          nullptr, SUNMatrixWrapper(x.getLength(), ubw, lbw, x.get_ctx())
+          nullptr, SUNMatrixWrapper(x.size(), ubw, lbw, x.get_ctx())
       ) {
     solver_
         = SUNLinSol_Band(const_cast<N_Vector>(x.getNVector()), A_, x.get_ctx());
@@ -191,7 +191,7 @@ SUNLinSolBand::SUNLinSolBand(AmiVector const& x, int ubw, int lbw)
 
 SUNLinSolDense::SUNLinSolDense(AmiVector const& x)
     : SUNLinSolWrapper(
-          nullptr, SUNMatrixWrapper(x.getLength(), x.getLength(), x.get_ctx())
+          nullptr, SUNMatrixWrapper(x.size(), x.size(), x.get_ctx())
       ) {
     solver_ = SUNLinSol_Dense(
         const_cast<N_Vector>(x.getNVector()), A_, x.get_ctx()
@@ -212,7 +212,7 @@ SUNLinSolKLU::SUNLinSolKLU(
     : SUNLinSolWrapper(
           nullptr,
           SUNMatrixWrapper(
-              x.getLength(), x.getLength(), nnz, sparsetype, x.get_ctx()
+              x.size(), x.size(), nnz, sparsetype, x.get_ctx()
           )
       ) {
     solver_
diff --git a/src/vector.cpp b/src/vector.cpp
index 07802c50db..2ccde593e1 100644
--- a/src/vector.cpp
+++ b/src/vector.cpp
@@ -21,7 +21,7 @@ const_N_Vector AmiVector::getNVector() const { return nvec_; }
 
 std::vector<realtype> const& AmiVector::getVector() const { return vec_; }
 
-int AmiVector::getLength() const { return gsl::narrow<int>(vec_.size()); }
+int AmiVector::size() const { return gsl::narrow<int>(vec_.size()); }
 
 void AmiVector::zero() { set(0.0); }
 
@@ -48,11 +48,11 @@ realtype const& AmiVector::at(int const pos) const {
 }
 
 void AmiVector::copy(AmiVector const& other) {
-    if (getLength() != other.getLength())
+    if (size() != other.size())
         throw AmiException(
             "Dimension of AmiVector (%i) does not "
             "match input dimension (%i)",
-            getLength(), other.getLength()
+            size(), other.size()
         );
     std::ranges::copy(other.vec_, vec_.begin());
 }
@@ -85,8 +85,8 @@ AmiVectorArray::AmiVectorArray(
 
 AmiVectorArray& AmiVectorArray::operator=(AmiVectorArray const& other) {
     vec_array_ = other.vec_array_;
-    nvec_array_.resize(other.getLength());
-    for (int idx = 0; idx < other.getLength(); idx++) {
+    nvec_array_.resize(other.size());
+    for (int idx = 0; idx < other.size(); idx++) {
         nvec_array_.at(idx) = vec_array_.at(idx).getNVector();
     }
     return *this;
@@ -94,8 +94,8 @@ AmiVectorArray& AmiVectorArray::operator=(AmiVectorArray const& other) {
 
 AmiVectorArray::AmiVectorArray(AmiVectorArray const& vaold)
     : vec_array_(vaold.vec_array_) {
-    nvec_array_.resize(vaold.getLength());
-    for (int idx = 0; idx < vaold.getLength(); idx++) {
+    nvec_array_.resize(vaold.size());
+    for (int idx = 0; idx < vaold.size(); idx++) {
         nvec_array_.at(idx) = vec_array_.at(idx).getNVector();
     }
 }
@@ -132,7 +132,7 @@ AmiVector const& AmiVectorArray::operator[](int const pos) const {
     return vec_array_.at(pos);
 }
 
-int AmiVectorArray::getLength() const {
+int AmiVectorArray::size() const {
     return gsl::narrow<int>(vec_array_.size());
 }
 
@@ -145,7 +145,7 @@ void AmiVectorArray::flatten_to_vector(std::vector<realtype>& vec) const {
     int n_outer = gsl::narrow<int>(vec_array_.size());
     if (n_outer == 0)
         return; // nothing to do ...
-    int n_inner = vec_array_.at(0).getLength();
+    int n_inner = vec_array_.at(0).size();
 
     if (gsl::narrow<int>(vec.size()) != n_inner * n_outer) {
         throw AmiException(
@@ -162,14 +162,14 @@ void AmiVectorArray::flatten_to_vector(std::vector<realtype>& vec) const {
 }
 
 void AmiVectorArray::copy(AmiVectorArray const& other) {
-    if (getLength() != other.getLength())
+    if (size() != other.size())
         throw AmiException(
             "Dimension of AmiVectorArray (%i) does not "
             "match input dimension (%i)",
-            getLength(), other.getLength()
+            size(), other.size()
         );
 
-    for (int iv = 0; iv < getLength(); ++iv) {
+    for (int iv = 0; iv < size(); ++iv) {
         vec_array_.at(iv).copy(other.vec_array_.at(iv));
         nvec_array_[iv] = vec_array_.at(iv).getNVector();
     }
diff --git a/tests/cpp/unittests/testMisc.cpp b/tests/cpp/unittests/testMisc.cpp
index 22d11d23cc..131c3c9ff5 100644
--- a/tests/cpp/unittests/testMisc.cpp
+++ b/tests/cpp/unittests/testMisc.cpp
@@ -472,7 +472,7 @@ TEST_F(AmiVectorTest, Vector) {
     AmiVector av2(nvec);
     ASSERT_NE(av.data(), av2.data());
 
-    for (int i = 0; i < av.getLength(); ++i) {
+    for (int i = 0; i < av.size(); ++i) {
         ASSERT_EQ(av.at(i), NV_Ith_S(nvec, i));
         ASSERT_EQ(av.at(i), av2.at(i));
     }
@@ -487,7 +487,7 @@ TEST_F(AmiVectorTest, VectorArray) {
     AmiVectorArray ava(4, 3, sunctx);
     AmiVector av1(vec1, sunctx), av2(vec2, sunctx), av3(vec3, sunctx);
     std::vector<AmiVector> avs{av1, av2, av3};
-    for (int i = 0; i < ava.getLength(); ++i)
+    for (int i = 0; i < ava.size(); ++i)
         ava[i] = avs.at(i);
 
     std::vector<double> badLengthVector(13, 0.0);
@@ -495,10 +495,10 @@ TEST_F(AmiVectorTest, VectorArray) {
 
     ASSERT_THROW(ava.flatten_to_vector(badLengthVector), AmiException);
     ava.flatten_to_vector(flattened);
-    for (int i = 0; i < ava.getLength(); ++i) {
+    for (int i = 0; i < ava.size(); ++i) {
         AmiVector const av = ava[i];
-        for (int j = 0; j < av.getLength(); ++j)
-            ASSERT_EQ(flattened.at(i * av.getLength() + j), av.at(j));
+        for (int j = 0; j < av.size(); ++j)
+            ASSERT_EQ(flattened.at(i * av.size() + j), av.at(j));
     }
 
     std::stringstream ss;