diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl
index 6f3840a..bbc9e49 100644
--- a/src/MOI_wrapper/MOI_wrapper.jl
+++ b/src/MOI_wrapper/MOI_wrapper.jl
@@ -2155,6 +2155,46 @@ function MOI.modify(
     return
 end
 
+function MOI.modify(
+    model::Optimizer,
+    cis::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S}},
+    changes::Vector{MOI.ScalarCoefficientChange{Float64}},
+) where {S}
+    nels = length(cis)
+    @assert nels == length(changes)
+    rows = Vector{Cint}(undef, nels)
+    cols = Vector{Cint}(undef, nels)
+    coefs = Vector{Cdouble}(undef, nels)
+    for i in 1:nels
+        rows[i] = Cint(_info(model, cis[i]).row)
+        cols[i] = column(model, changes[i].variable)
+        coefs[i] = changes[i].new_coefficient
+    end
+    for row in unique(rows)
+        nnz = glp_get_mat_row(model, row, C_NULL, C_NULL)
+        indices, coefficients = zeros(Cint, nnz), zeros(Cdouble, nnz)
+        glp_get_mat_row(model, row, offset(indices), offset(coefficients))
+        idxs_changed_in_row = findall(x -> x == row, rows)
+        for i in idxs_changed_in_row
+            index = something(findfirst(isequal(cols[i]), indices), 0)
+            if index > 0
+                coefficients[index] = coefs[i]
+            else
+                push!(indices, cols[i])
+                push!(coefficients, coefs[i])
+            end
+        end
+        glp_set_mat_row(
+            model,
+            row,
+            length(indices),
+            offset(indices),
+            offset(coefficients),
+        )
+    end
+    return
+end
+
 function MOI.set(
     model::Optimizer,
     ::MOI.ConstraintFunction,
diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl
index be5a7f1..3fe5bd0 100644
--- a/test/MOI_wrapper.jl
+++ b/test/MOI_wrapper.jl
@@ -544,6 +544,51 @@ function test_variable_basis_status()
     return
 end
 
+function test_multiple_modifications()
+    model = GLPK.Optimizer()
+
+    x = MOI.add_variables(model, 3)
+
+    saf = MOI.ScalarAffineFunction(
+        [
+            MOI.ScalarAffineTerm(1.0, x[1]),
+            MOI.ScalarAffineTerm(1.0, x[2]),
+            MOI.ScalarAffineTerm(1.0, x[3]),
+        ],
+        0.0,
+    )
+    ci1 = MOI.add_constraint(model, saf, MOI.LessThan(1.0))
+    ci2 = MOI.add_constraint(model, saf, MOI.LessThan(2.0))
+
+    MOI.set(
+        model,
+        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
+        saf,
+    )
+
+    fc1 = MOI.get(model, MOI.ConstraintFunction(), ci1)
+    @test MOI.coefficient.(fc1.terms) == [1.0, 1.0, 1.0]
+    fc2 = MOI.get(model, MOI.ConstraintFunction(), ci2)
+    @test MOI.coefficient.(fc2.terms) == [1.0, 1.0, 1.0]
+    obj = MOI.get(
+        model,
+        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
+    )
+    @test MOI.coefficient.(obj.terms) == [1.0, 1.0, 1.0]
+
+    changes_cis = [
+        MOI.ScalarCoefficientChange(MOI.VariableIndex(1), 4.0)
+        MOI.ScalarCoefficientChange(MOI.VariableIndex(1), 0.5)
+        MOI.ScalarCoefficientChange(MOI.VariableIndex(3), 2.0)
+    ]
+    MOI.modify(model, [ci1, ci2, ci2], changes_cis)
+
+    fc1 = MOI.get(model, MOI.ConstraintFunction(), ci1)
+    @test MOI.coefficient.(fc1.terms) == [4.0, 1.0, 1.0]
+    fc2 = MOI.get(model, MOI.ConstraintFunction(), ci2)
+    @test MOI.coefficient.(fc2.terms) == [0.5, 1.0, 2.0]
+end
+
 end  # module
 
 TestMOIWrapper.runtests()