diff --git a/src/nlp_expr.jl b/src/nlp_expr.jl
index 67354a1cfb7..61e018588d0 100644
--- a/src/nlp_expr.jl
+++ b/src/nlp_expr.jl
@@ -634,6 +634,33 @@ function moi_function(model::GenericModel, f::GenericNonlinearExpr{V}) where {V}
     return ret
 end
 
+# A backwards-compatible function to preserve behavior prior to #4032. This was
+# used by Plasmo.jl
+function moi_function(f::GenericNonlinearExpr{V}) where {V}
+    ret = MOI.ScalarNonlinearFunction(f.head, similar(f.args))
+    stack = Tuple{MOI.ScalarNonlinearFunction,Int,GenericNonlinearExpr{V}}[]
+    for i in length(f.args):-1:1
+        if f.args[i] isa GenericNonlinearExpr{V}
+            push!(stack, (ret, i, f.args[i]))
+        else
+            ret.args[i] = moi_function(f.args[i])
+        end
+    end
+    while !isempty(stack)
+        parent, i, arg = pop!(stack)
+        child = MOI.ScalarNonlinearFunction(arg.head, similar(arg.args))
+        parent.args[i] = child
+        for j in length(arg.args):-1:1
+            if arg.args[j] isa GenericNonlinearExpr{V}
+                push!(stack, (child, j, arg.args[j]))
+            else
+                child.args[j] = moi_function(arg.args[j])
+            end
+        end
+    end
+    return ret
+end
+
 jump_function(::GenericModel{T}, x::Number) where {T} = convert(T, x)
 
 function jump_function(model::GenericModel, f::MOI.ScalarNonlinearFunction)
diff --git a/test/test_nlp_expr.jl b/test/test_nlp_expr.jl
index bbfe010ebf2..84859bea94d 100644
--- a/test/test_nlp_expr.jl
+++ b/test/test_nlp_expr.jl
@@ -1381,4 +1381,37 @@ function test_custom_array()
     return
 end
 
+function test_scalar_nonlinear_moi_function()
+    model = Model()
+    @variable(model, x)
+    y1 = atan(x, 2)
+    y1_moi = MOI.ScalarNonlinearFunction(:atan, Any[index(x), 2])
+    y2 = y1 + y1
+    y2_moi = MOI.ScalarNonlinearFunction(:+, Any[y1_moi, y1_moi])
+    y3 = exp(y2)
+    y3_moi = MOI.ScalarNonlinearFunction(:exp, Any[y2_moi])
+    # Test y1
+    @test isapprox(moi_function(y1), y1_moi)
+    @test isempty(model.subexpressions)
+    @test isapprox(moi_function(model, y1), y1_moi)
+    @test length(model.subexpressions) == 1
+    @test isapprox(model.subexpressions[objectid(y1)], y1_moi)
+    # Test y2
+    @test isapprox(moi_function(y2), y2_moi)
+    @test length(model.subexpressions) == 1
+    @test isapprox(moi_function(model, y2), y2_moi)
+    @test length(model.subexpressions) == 2
+    @test isapprox(model.subexpressions[objectid(y1)], y1_moi)
+    @test isapprox(model.subexpressions[objectid(y2)], y2_moi)
+    # Test y3
+    @test isapprox(moi_function(y3), y3_moi)
+    @test length(model.subexpressions) == 2
+    @test isapprox(moi_function(model, y3), y3_moi)
+    @test length(model.subexpressions) == 3
+    @test isapprox(model.subexpressions[objectid(y1)], y1_moi)
+    @test isapprox(model.subexpressions[objectid(y2)], y2_moi)
+    @test isapprox(model.subexpressions[objectid(y3)], y3_moi)
+    return
+end
+
 end  # module