diff --git a/src/Utilities/matrix_of_constraints.jl b/src/Utilities/matrix_of_constraints.jl
index 43d4f80448..69e6ec9683 100644
--- a/src/Utilities/matrix_of_constraints.jl
+++ b/src/Utilities/matrix_of_constraints.jl
@@ -739,6 +739,32 @@ function MOI.modify(
     return
 end
 
+function MOI.modify(
+    model::MatrixOfConstraints,
+    ci::MOI.ConstraintIndex,
+    change::MOI.MultirowChange,
+)
+    r = rows(model, ci)
+    for (output_index, new_coefficient) in change.new_coefficients
+        if !modify_coefficients(
+            model.coefficients,
+            r[output_index],
+            change.variable.value,
+            new_coefficient,
+        )
+            throw(
+                MOI.ModifyConstraintNotAllowed(
+                    ci,
+                    change,
+                    "cannot set a new non-zero coefficient because no entry " *
+                    "exists in the sparse matrix of `MatrixOfConstraints`",
+                ),
+            )
+        end
+    end
+    return
+end
+
 # See https://github.com/jump-dev/MathOptInterface.jl/pull/2976
 # Ideally we would have made it so that `modify_constants` operated like
 # `modify_coefficients` and returned a `Bool` indicating success. But we didn't,
diff --git a/src/Utilities/sparse_matrix.jl b/src/Utilities/sparse_matrix.jl
index 45b831c03b..cc9ad4a287 100644
--- a/src/Utilities/sparse_matrix.jl
+++ b/src/Utilities/sparse_matrix.jl
@@ -239,7 +239,8 @@ function extract_function(
             continue
         end
         r = _shift(A.rowval[idx], _indexing(A), OneBasedIndexing())
-        if r == row
+        # `modify_coefficients` can create zeros
+        if r == row && !iszero(A.nzval[idx])
             push!(
                 func.terms,
                 MOI.ScalarAffineTerm(A.nzval[idx], MOI.VariableIndex(col)),
@@ -266,10 +267,11 @@ function _extract_column_as_function(
     func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{T}[], constant)
     for i in SparseArrays.nzrange(A, col)
         row = _shift(A.rowval[i], _indexing(A), OneBasedIndexing())
-        push!(
-            func.terms,
-            MOI.ScalarAffineTerm(value_map(A.nzval[i]), MOI.VariableIndex(row)),
-        )
+        val = value_map(A.nzval[i])
+        # `modify_coefficients` can create zeros
+        if !iszero(val)
+            push!(func.terms, MOI.ScalarAffineTerm(val, MOI.VariableIndex(row)))
+        end
     end
     return func
 end
@@ -293,16 +295,19 @@ function extract_function(
             if row != rows[output_index]
                 continue
             end
-            push!(
-                func.terms,
-                MOI.VectorAffineTerm(
-                    output_index,
-                    MOI.ScalarAffineTerm(
-                        A.nzval[idx[col]],
-                        MOI.VariableIndex(col),
+            # `modify_coefficients` can create zeros
+            if !iszero(A.nzval[idx[col]])
+                push!(
+                    func.terms,
+                    MOI.VectorAffineTerm(
+                        output_index,
+                        MOI.ScalarAffineTerm(
+                            A.nzval[idx[col]],
+                            MOI.VariableIndex(col),
+                        ),
                     ),
-                ),
-            )
+                )
+            end
             idx[col] += 1
         end
     end
diff --git a/test/Utilities/test_matrix_of_constraints.jl b/test/Utilities/test_matrix_of_constraints.jl
index f545c2b3b9..c6be87a5fa 100644
--- a/test/Utilities/test_matrix_of_constraints.jl
+++ b/test/Utilities/test_matrix_of_constraints.jl
@@ -665,6 +665,7 @@ function test_modify_scalar_coefficient_change()
     @test f ≈ 5x + 3y
     MOI.modify(model, c, MOI.ScalarCoefficientChange(y, 0))
     f = MOI.get(model, MOI.ConstraintFunction(), c)
+    @test MOI.Utilities.is_canonical(f)
     @test f ≈ 5x + 0y
     return
 end
@@ -826,6 +827,124 @@ function test_modify_set_constants()
     return
 end
 
+function test_modify_multirow_change()
+    model = _new_VectorSets()
+    src = MOI.Utilities.Model{Int}()
+    x = MOI.add_variables(src, 2)
+    terms = [
+        MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2, x[1])),
+        MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(3, x[2])),
+        MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(4, x[1])),
+        MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(5, x[2])),
+    ]
+    func = MOI.VectorAffineFunction(terms, [0, 0])
+    c = MOI.add_constraint(src, func, MOI.Nonnegatives(2))
+    index_map = MOI.copy_to(model, src)
+    ci = index_map[c]
+    f = MOI.get(model, MOI.ConstraintFunction(), ci)
+    @test length(f.terms) == 4
+    x1 = index_map[x[1]]
+    x2 = index_map[x[2]]
+    MOI.modify(model, ci, MOI.MultirowChange(x1, [(1, 7), (2, 8)]))
+    f = MOI.get(model, MOI.ConstraintFunction(), ci)
+    coefs = Dict(
+        (t.output_index, t.scalar_term.variable) =>
+            t.scalar_term.coefficient for t in f.terms
+    )
+    @test coefs[(1, x1)] == 7
+    @test coefs[(2, x1)] == 8
+    @test coefs[(1, x2)] == 3
+    @test coefs[(2, x2)] == 5
+    return
+end
+
+function test_modify_multirow_change_single_row()
+    model = _new_VectorSets()
+    src = MOI.Utilities.Model{Int}()
+    x = MOI.add_variables(src, 2)
+    terms = [
+        MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2, x[1])),
+        MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(3, x[2])),
+        MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(4, x[1])),
+        MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(5, x[2])),
+    ]
+    func = MOI.VectorAffineFunction(terms, [0, 0])
+    c = MOI.add_constraint(src, func, MOI.Nonnegatives(2))
+    index_map = MOI.copy_to(model, src)
+    ci = index_map[c]
+    x2 = index_map[x[2]]
+    MOI.modify(model, ci, MOI.MultirowChange(x2, [(2, 9)]))
+    f = MOI.get(model, MOI.ConstraintFunction(), ci)
+    coefs = Dict(
+        (t.output_index, t.scalar_term.variable) =>
+            t.scalar_term.coefficient for t in f.terms
+    )
+    x1 = index_map[x[1]]
+    @test coefs[(1, x1)] == 2
+    @test coefs[(2, x1)] == 4
+    @test coefs[(1, x2)] == 3
+    @test coefs[(2, x2)] == 9
+    return
+end
+
+function test_modify_multirow_change_to_zero()
+    model = _new_VectorSets()
+    src = MOI.Utilities.Model{Int}()
+    x = MOI.add_variables(src, 2)
+    terms = [
+        MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2, x[1])),
+        MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(3, x[2])),
+        MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(4, x[1])),
+        MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(5, x[2])),
+    ]
+    func = MOI.VectorAffineFunction(terms, [0, 0])
+    c = MOI.add_constraint(src, func, MOI.Nonnegatives(2))
+    index_map = MOI.copy_to(model, src)
+    ci = index_map[c]
+    x1 = index_map[x[1]]
+    x2 = index_map[x[2]]
+    MOI.modify(model, ci, MOI.MultirowChange(x1, [(1, 0)]))
+    f = MOI.get(model, MOI.ConstraintFunction(), ci)
+    @test MOI.Utilities.is_canonical(f)
+    @test length(f.terms) == 3
+    coefs = Dict(
+        (t.output_index, t.scalar_term.variable) =>
+            t.scalar_term.coefficient for t in f.terms
+    )
+    @test !haskey(coefs, (1, x1))
+    @test coefs[(2, x1)] == 4
+    @test coefs[(1, x2)] == 3
+    @test coefs[(2, x2)] == 5
+    return
+end
+
+function test_modify_multirow_change_no_entry()
+    model = _new_VectorSets()
+    src = MOI.Utilities.Model{Int}()
+    x = MOI.add_variables(src, 3)
+    terms = [
+        MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2, x[1])),
+        MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(4, x[1])),
+    ]
+    func = MOI.VectorAffineFunction(terms, [0, 0])
+    c = MOI.add_constraint(src, func, MOI.Nonnegatives(2))
+    index_map = MOI.copy_to(model, src)
+    ci = index_map[c]
+    x3 = index_map[x[3]]
+    MOI.modify(model, ci, MOI.MultirowChange(x3, [(1, 0)]))
+    change = MOI.MultirowChange(x3, [(1, 5)])
+    @test_throws(
+        MOI.ModifyConstraintNotAllowed(
+            ci,
+            change,
+            "cannot set a new non-zero coefficient because no entry " *
+            "exists in the sparse matrix of `MatrixOfConstraints`",
+        ),
+        MOI.modify(model, ci, change),
+    )
+    return
+end
+
 function test_unsupported_constraint()
     model = _new_ScalarSets()
     x = MOI.VariableIndex(1)