diff --git a/src/parameters.jl b/src/parameters.jl
index c7a484b1..d8e36962 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -266,6 +266,222 @@ function _quadratic_objective_set_forward!(model::POI.Optimizer{T}) where {T}
     return
 end
 
+function _cubic_objective_set_forward!(model::POI.Optimizer{T}) where {T}
+    pf = MOI.get(
+        model,
+        POI.ParametricObjectiveFunction{POI.ParametricCubicFunction{T}}(),
+    )
+    sensitivity_data = _get_sensitivity_data(model)
+
+    cte = zero(T)
+    affine_terms = MOI.ScalarAffineTerm{T}[]
+    quadratic_terms = MOI.ScalarQuadraticTerm{T}[]
+
+    # pvv terms: Δp * coeff → quadratic perturbation
+    for term in POI._cubic_pvv_terms(pf)
+        p = term.index_1
+        v1 = term.index_2
+        v2 = term.index_3
+        Δp = get(sensitivity_data.parameter_input_forward, p, zero(T))
+        if !iszero(Δp)
+            qcoeff = Δp * term.coefficient
+            if v1 == v2
+                qcoeff *= 2  # MOI diagonal convention
+            end
+            push!(
+                quadratic_terms,
+                MOI.ScalarQuadraticTerm{T}(qcoeff, v1, v2),
+            )
+        end
+    end
+
+    # ppv terms: chain rule on two params → affine perturbation
+    for term in POI._cubic_ppv_terms(pf)
+        p1 = term.index_1
+        p2 = term.index_2
+        v = term.index_3
+        Δp1 = get(sensitivity_data.parameter_input_forward, p1, zero(T))
+        Δp2 = get(sensitivity_data.parameter_input_forward, p2, zero(T))
+        p1_val = MOI.get(model, MOI.VariablePrimal(), p1)
+        p2_val = MOI.get(model, MOI.VariablePrimal(), p2)
+        acoeff =
+            Δp1 * term.coefficient * p2_val + Δp2 * term.coefficient * p1_val
+        if !iszero(acoeff)
+            push!(affine_terms, MOI.ScalarAffineTerm{T}(acoeff, v))
+        end
+    end
+
+    # ppp terms: chain rule on three params → constant perturbation
+    for term in POI._cubic_ppp_terms(pf)
+        p1 = term.index_1
+        p2 = term.index_2
+        p3 = term.index_3
+        Δp1 = get(sensitivity_data.parameter_input_forward, p1, zero(T))
+        Δp2 = get(sensitivity_data.parameter_input_forward, p2, zero(T))
+        Δp3 = get(sensitivity_data.parameter_input_forward, p3, zero(T))
+        p1_val = MOI.get(model, MOI.VariablePrimal(), p1)
+        p2_val = MOI.get(model, MOI.VariablePrimal(), p2)
+        p3_val = MOI.get(model, MOI.VariablePrimal(), p3)
+        cte += term.coefficient * (
+            Δp1 * p2_val * p3_val +
+            p1_val * Δp2 * p3_val +
+            p1_val * p2_val * Δp3
+        )
+    end
+
+    # Degree-2: p terms (affine parameter → constant perturbation)
+    for term in pf.p
+        p = term.variable
+        Δp = get(sensitivity_data.parameter_input_forward, p, zero(T))
+        cte += Δp * term.coefficient
+    end
+    # Degree-2: pp terms (parameter-parameter → constant perturbation)
+    for term in pf.pp
+        p_1 = term.variable_1
+        p_2 = term.variable_2
+        Δp1 = get(sensitivity_data.parameter_input_forward, p_1, zero(T))
+        Δp2 = get(sensitivity_data.parameter_input_forward, p_2, zero(T))
+        p1_val = MOI.get(model, MOI.VariablePrimal(), p_1)
+        p2_val = MOI.get(model, MOI.VariablePrimal(), p_2)
+        cte +=
+            Δp1 *
+            term.coefficient *
+            p2_val /
+            ifelse(p_1 === p_2, T(2), T(1))
+        cte +=
+            Δp2 *
+            term.coefficient *
+            p1_val /
+            ifelse(p_1 === p_2, T(2), T(1))
+    end
+    # Degree-2: pv terms (parameter-variable → affine perturbation)
+    for term in pf.pv
+        p = term.variable_1
+        Δp = get(sensitivity_data.parameter_input_forward, p, zero(T))
+        if !iszero(Δp)
+            push!(
+                affine_terms,
+                MOI.ScalarAffineTerm{T}(
+                    Δp * term.coefficient,
+                    term.variable_2,
+                ),
+            )
+        end
+    end
+
+    # Send perturbation to inner model
+    if !isempty(quadratic_terms) || !isempty(affine_terms) || !iszero(cte)
+        MOI.set(
+            model.optimizer,
+            ForwardObjectiveFunction(),
+            MOI.ScalarQuadraticFunction{T}(quadratic_terms, affine_terms, cte),
+        )
+    end
+    return
+end
+
+function _cubic_objective_get_reverse!(model::POI.Optimizer{T}) where {T}
+    pf = MOI.get(
+        model,
+        POI.ParametricObjectiveFunction{POI.ParametricCubicFunction{T}}(),
+    )
+    pvv_terms = POI._cubic_pvv_terms(pf)
+    ppv_terms = POI._cubic_ppv_terms(pf)
+    ppp_terms = POI._cubic_ppp_terms(pf)
+    p_terms = pf.p
+    pp_terms = pf.pp
+    pv_terms = pf.pv
+    if isempty(pvv_terms) && isempty(ppv_terms) && isempty(ppp_terms) &&
+       isempty(p_terms) && isempty(pp_terms) && isempty(pv_terms)
+        return
+    end
+    sensitivity_data = _get_sensitivity_data(model)
+    grad_pf = MOI.get(model.optimizer, ReverseObjectiveFunction())
+    d_cte = MOI.constant(grad_pf)
+
+    # pvv terms: ∇p += 2 * coeff * dQ[v1, v2]
+    for term in pvv_terms
+        p = term.index_1
+        v1 = term.index_2
+        v2 = term.index_3
+        dQ_ij = quad_sym_half(grad_pf, v1, v2)
+        value = get!(sensitivity_data.parameter_output_backward, p, zero(T))
+        sensitivity_data.parameter_output_backward[p] =
+            value + T(2) * term.coefficient * dQ_ij
+    end
+
+    # ppv terms: chain rule on two params
+    for term in ppv_terms
+        p1 = term.index_1
+        p2 = term.index_2
+        v = term.index_3
+        dq_v = JuMP.coefficient(grad_pf, v)
+        p1_val = MOI.get(model, MOI.VariablePrimal(), p1)
+        p2_val = MOI.get(model, MOI.VariablePrimal(), p2)
+        val1 = get!(sensitivity_data.parameter_output_backward, p1, zero(T))
+        val2 = get!(sensitivity_data.parameter_output_backward, p2, zero(T))
+        sensitivity_data.parameter_output_backward[p1] =
+            val1 + term.coefficient * p2_val * dq_v
+        sensitivity_data.parameter_output_backward[p2] =
+            val2 + term.coefficient * p1_val * dq_v
+    end
+
+    # ppp terms: chain rule on three params
+    for term in ppp_terms
+        p1 = term.index_1
+        p2 = term.index_2
+        p3 = term.index_3
+        p1_val = MOI.get(model, MOI.VariablePrimal(), p1)
+        p2_val = MOI.get(model, MOI.VariablePrimal(), p2)
+        p3_val = MOI.get(model, MOI.VariablePrimal(), p3)
+        val1 = get!(sensitivity_data.parameter_output_backward, p1, zero(T))
+        val2 = get!(sensitivity_data.parameter_output_backward, p2, zero(T))
+        val3 = get!(sensitivity_data.parameter_output_backward, p3, zero(T))
+        sensitivity_data.parameter_output_backward[p1] =
+            val1 + term.coefficient * p2_val * p3_val * d_cte
+        sensitivity_data.parameter_output_backward[p2] =
+            val2 + term.coefficient * p1_val * p3_val * d_cte
+        sensitivity_data.parameter_output_backward[p3] =
+            val3 + term.coefficient * p1_val * p2_val * d_cte
+    end
+
+    # Degree-2: p terms
+    for term in p_terms
+        p = term.variable
+        value = get!(sensitivity_data.parameter_output_backward, p, zero(T))
+        sensitivity_data.parameter_output_backward[p] =
+            value + term.coefficient * d_cte
+    end
+    # Degree-2: pp terms
+    for term in pp_terms
+        p_1 = term.variable_1
+        p_2 = term.variable_2
+        value_1 = get!(sensitivity_data.parameter_output_backward, p_1, zero(T))
+        value_2 = get!(sensitivity_data.parameter_output_backward, p_2, zero(T))
+        sensitivity_data.parameter_output_backward[p_1] =
+            value_1 +
+            term.coefficient *
+            d_cte *
+            MOI.get(model, MOI.VariablePrimal(), p_2) /
+            ifelse(p_1 === p_2, T(2), T(1))
+        sensitivity_data.parameter_output_backward[p_2] =
+            value_2 +
+            term.coefficient *
+            d_cte *
+            MOI.get(model, MOI.VariablePrimal(), p_1) /
+            ifelse(p_1 === p_2, T(2), T(1))
+    end
+    # Degree-2: pv terms
+    for term in pv_terms
+        p = term.variable_1
+        v = term.variable_2
+        value = get!(sensitivity_data.parameter_output_backward, p, zero(T))
+        sensitivity_data.parameter_output_backward[p] =
+            value + term.coefficient * JuMP.coefficient(grad_pf, v)
+    end
+    return
+end
+
 function empty_input_sensitivities!(model::POI.Optimizer{T}) where {T}
     empty_input_sensitivities!(model.optimizer)
     model.ext[_SENSITIVITY_DATA] = SensitivityData{T}()
@@ -287,6 +503,8 @@ function forward_differentiate!(model::POI.Optimizer{T}) where {T}
         _affine_objective_set_forward!(model)
     elseif obj_type <: POI.ParametricQuadraticFunction
         _quadratic_objective_set_forward!(model)
+    elseif obj_type <: POI.ParametricCubicFunction
+        _cubic_objective_set_forward!(model)
     end
     forward_differentiate!(model.optimizer)
     return
@@ -511,6 +729,8 @@ function reverse_differentiate!(model::POI.Optimizer)
         _affine_objective_get_reverse!(model)
     elseif obj_type <: POI.ParametricQuadraticFunction
         _quadratic_objective_get_reverse!(model)
+    elseif obj_type <: POI.ParametricCubicFunction
+        _cubic_objective_get_reverse!(model)
     end
     return
 end
diff --git a/test/parameters.jl b/test/parameters.jl
index 9b4f1968..58ade75d 100644
--- a/test/parameters.jl
+++ b/test/parameters.jl
@@ -710,6 +710,379 @@ function test_empty_cache()
     return
 end
 
+# ============================================================
+# Tests for cubic objective differentiation (pvv, ppv, ppp)
+# ============================================================
+
+function test_cubic_obj_pvv_nonbinding()
+    # min p * x^2 + x  s.t. x >= -100
+    # FOC: 2p*x + 1 = 0  → x = -1/(2p)
+    # dx/dp = 1/(2p^2)
+    model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer))
+    set_silent(model)
+    @variable(model, x)
+    @variable(model, p in Parameter(2.0))
+    @constraint(model, x >= -100)
+    @objective(model, Min, p * x * x + x)
+    for p_val in [1.0, 2.0, 3.0]
+        set_parameter_value(p, p_val)
+        optimize!(model)
+        @test isapprox(value(x), -1 / (2 * p_val), atol = 1e-6)
+        for dir_p in [1.0, 2.0]
+            DiffOpt.set_forward_parameter(model, p, dir_p)
+            DiffOpt.forward_differentiate!(model)
+            @test isapprox(
+                DiffOpt.get_forward_variable(model, x),
+                dir_p / (2 * p_val^2),
+                atol = 1e-5,
+            )
+        end
+        for dir_x in [1.0, 2.0]
+            DiffOpt.set_reverse_variable(model, x, dir_x)
+            DiffOpt.reverse_differentiate!(model)
+            @test isapprox(
+                DiffOpt.get_reverse_parameter(model, p),
+                dir_x / (2 * p_val^2),
+                atol = 1e-5,
+            )
+        end
+    end
+    return
+end
+
+function test_cubic_obj_pvv_two_vars()
+    # min p*x*y + x^2 + y^2 + x + y  s.t. x >= -10, y >= -10
+    # FOC: p*y + 2x + 1 = 0, p*x + 2y + 1 = 0
+    # By symmetry x = y, so: p*x + 2x + 1 = 0 → x = -1/(2+p)
+    # dx/dp = 1/(2+p)^2
+    model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer))
+    set_silent(model)
+    @variable(model, x)
+    @variable(model, y)
+    @variable(model, p in Parameter(1.0))
+    @constraint(model, x >= -10)
+    @constraint(model, y >= -10)
+    @objective(model, Min, p * x * y + x^2 + y^2 + x + y)
+    for p_val in [0.5, 1.0, 1.5]
+        set_parameter_value(p, p_val)
+        optimize!(model)
+        expected_xy = -1 / (2 + p_val)
+        @test isapprox(value(x), expected_xy, atol = 1e-5)
+        @test isapprox(value(y), expected_xy, atol = 1e-5)
+        for dir_p in [1.0, 2.0]
+            DiffOpt.set_forward_parameter(model, p, dir_p)
+            DiffOpt.forward_differentiate!(model)
+            dx_dp = dir_p / (2 + p_val)^2
+            @test isapprox(
+                DiffOpt.get_forward_variable(model, x),
+                dx_dp,
+                atol = 1e-5,
+            )
+            @test isapprox(
+                DiffOpt.get_forward_variable(model, y),
+                dx_dp,
+                atol = 1e-5,
+            )
+        end
+        for dir_x in [1.0, 2.0]
+            DiffOpt.set_reverse_variable(model, x, dir_x)
+            DiffOpt.set_reverse_variable(model, y, 0.0)
+            DiffOpt.reverse_differentiate!(model)
+            @test isapprox(
+                DiffOpt.get_reverse_parameter(model, p),
+                dir_x / (2 + p_val)^2,
+                atol = 1e-5,
+            )
+        end
+    end
+    return
+end
+
+function test_cubic_obj_pvv_binding()
+    # min p * x^2  s.t. x >= 1
+    # x* = 1 (binding), dx/dp = 0
+    model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer))
+    set_silent(model)
+    @variable(model, x)
+    @variable(model, p in Parameter(2.0))
+    @constraint(model, x >= 1)
+    @objective(model, Min, p * x * x)
+    for p_val in [1.0, 2.0, 3.0]
+        set_parameter_value(p, p_val)
+        optimize!(model)
+        @test isapprox(value(x), 1.0, atol = 1e-6)
+        DiffOpt.set_forward_parameter(model, p, 1.0)
+        DiffOpt.forward_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_forward_variable(model, x),
+            0.0,
+            atol = 1e-5,
+        )
+        DiffOpt.set_reverse_variable(model, x, 1.0)
+        DiffOpt.reverse_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_reverse_parameter(model, p),
+            0.0,
+            atol = 1e-5,
+        )
+    end
+    return
+end
+
+function test_cubic_obj_ppv()
+    # min p * q * x + x^2  s.t. x >= -10
+    # FOC: p*q + 2x = 0 → x = -p*q/2
+    # dx/dp = -q/2, dx/dq = -p/2
+    model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer))
+    set_silent(model)
+    @variable(model, x)
+    @variable(model, p in Parameter(2.0))
+    @variable(model, q in Parameter(3.0))
+    @constraint(model, x >= -100)
+    @objective(model, Min, p * q * x + x^2)
+    for p_val in [1.0, 2.0], q_val in [1.0, 3.0]
+        set_parameter_value(p, p_val)
+        set_parameter_value(q, q_val)
+        optimize!(model)
+        @test isapprox(value(x), -p_val * q_val / 2, atol = 1e-5)
+        # Forward dx/dp
+        DiffOpt.set_forward_parameter(model, p, 1.0)
+        DiffOpt.set_forward_parameter(model, q, 0.0)
+        DiffOpt.forward_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_forward_variable(model, x),
+            -q_val / 2,
+            atol = 1e-5,
+        )
+        # Forward dx/dq
+        DiffOpt.set_forward_parameter(model, p, 0.0)
+        DiffOpt.set_forward_parameter(model, q, 1.0)
+        DiffOpt.forward_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_forward_variable(model, x),
+            -p_val / 2,
+            atol = 1e-5,
+        )
+        # Reverse
+        DiffOpt.set_reverse_variable(model, x, 1.0)
+        DiffOpt.reverse_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_reverse_parameter(model, p),
+            -q_val / 2,
+            atol = 1e-5,
+        )
+        @test isapprox(
+            DiffOpt.get_reverse_parameter(model, q),
+            -p_val / 2,
+            atol = 1e-5,
+        )
+    end
+    return
+end
+
+function test_cubic_obj_mixed()
+    # min p * x^2 + q * x + r  s.t. x >= -100
+    # FOC: 2p*x + q = 0 → x = -q/(2p)
+    # dx/dp = q/(2p^2), dx/dq = -1/(2p), dx/dr = 0
+    model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer))
+    set_silent(model)
+    @variable(model, x)
+    @variable(model, p in Parameter(2.0))
+    @variable(model, q in Parameter(1.0))
+    @variable(model, r in Parameter(5.0))
+    @constraint(model, x >= -100)
+    @objective(model, Min, p * x * x + q * x + r)
+    for p_val in [1.0, 2.0], q_val in [1.0, 3.0], r_val in [1.0, 5.0]
+        set_parameter_value(p, p_val)
+        set_parameter_value(q, q_val)
+        set_parameter_value(r, r_val)
+        optimize!(model)
+        @test isapprox(value(x), -q_val / (2 * p_val), atol = 1e-5)
+        # dx/dp
+        DiffOpt.set_forward_parameter(model, p, 1.0)
+        DiffOpt.set_forward_parameter(model, q, 0.0)
+        DiffOpt.set_forward_parameter(model, r, 0.0)
+        DiffOpt.forward_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_forward_variable(model, x),
+            q_val / (2 * p_val^2),
+            atol = 1e-5,
+        )
+        # dx/dq
+        DiffOpt.set_forward_parameter(model, p, 0.0)
+        DiffOpt.set_forward_parameter(model, q, 1.0)
+        DiffOpt.forward_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_forward_variable(model, x),
+            -1 / (2 * p_val),
+            atol = 1e-5,
+        )
+        # dx/dr = 0
+        DiffOpt.set_forward_parameter(model, q, 0.0)
+        DiffOpt.set_forward_parameter(model, r, 1.0)
+        DiffOpt.forward_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_forward_variable(model, x),
+            0.0,
+            atol = 1e-5,
+        )
+        # Reverse
+        DiffOpt.set_reverse_variable(model, x, 1.0)
+        DiffOpt.reverse_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_reverse_parameter(model, p),
+            q_val / (2 * p_val^2),
+            atol = 1e-5,
+        )
+        @test isapprox(
+            DiffOpt.get_reverse_parameter(model, q),
+            -1 / (2 * p_val),
+            atol = 1e-5,
+        )
+        @test isapprox(
+            DiffOpt.get_reverse_parameter(model, r),
+            0.0,
+            atol = 1e-5,
+        )
+    end
+    return
+end
+
+function test_cubic_obj_multi_pvv()
+    # min p*x^2 + q*y^2 + x + y  s.t. x >= -10, y >= -10
+    # FOC: 2p*x + 1 = 0, 2q*y + 1 = 0
+    # x = -1/(2p), y = -1/(2q)
+    # dx/dp = 1/(2p^2), dx/dq = 0, dy/dp = 0, dy/dq = 1/(2q^2)
+    model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer))
+    set_silent(model)
+    @variable(model, x)
+    @variable(model, y)
+    @variable(model, p in Parameter(2.0))
+    @variable(model, q in Parameter(3.0))
+    @constraint(model, x >= -10)
+    @constraint(model, y >= -10)
+    @objective(model, Min, p * x * x + q * y * y + x + y)
+    for p_val in [1.0, 2.0], q_val in [2.0, 3.0]
+        set_parameter_value(p, p_val)
+        set_parameter_value(q, q_val)
+        optimize!(model)
+        @test isapprox(value(x), -1 / (2 * p_val), atol = 1e-5)
+        @test isapprox(value(y), -1 / (2 * q_val), atol = 1e-5)
+        # dx/dp
+        DiffOpt.set_forward_parameter(model, p, 1.0)
+        DiffOpt.set_forward_parameter(model, q, 0.0)
+        DiffOpt.forward_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_forward_variable(model, x),
+            1 / (2 * p_val^2),
+            atol = 1e-5,
+        )
+        @test isapprox(
+            DiffOpt.get_forward_variable(model, y),
+            0.0,
+            atol = 1e-5,
+        )
+        # Reverse: ∇x=1, ∇y=0
+        DiffOpt.set_reverse_variable(model, x, 1.0)
+        DiffOpt.set_reverse_variable(model, y, 0.0)
+        DiffOpt.reverse_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_reverse_parameter(model, p),
+            1 / (2 * p_val^2),
+            atol = 1e-5,
+        )
+        @test isapprox(
+            DiffOpt.get_reverse_parameter(model, q),
+            0.0,
+            atol = 1e-5,
+        )
+    end
+    return
+end
+
+function test_cubic_obj_ppv_and_pv()
+    # min p * q * x + r * x + x^2  s.t. x >= -100
+    # After substitution: (p*q + r) * x + x^2
+    # FOC: p*q + r + 2x = 0  → x = -(p*q + r)/2
+    # dx/dp = -q/2, dx/dq = -p/2, dx/dr = -1/2
+    model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer))
+    set_silent(model)
+    @variable(model, x)
+    @variable(model, p in Parameter(2.0))
+    @variable(model, q in Parameter(3.0))
+    @variable(model, r in Parameter(1.0))
+    @constraint(model, x >= -100)
+    @objective(model, Min, p * q * x + r * x + x^2)
+    for p_val in [1.0, 2.0], q_val in [1.0, 3.0], r_val in [0.5, 2.0]
+        set_parameter_value(p, p_val)
+        set_parameter_value(q, q_val)
+        set_parameter_value(r, r_val)
+        optimize!(model)
+        @test isapprox(
+            value(x),
+            -(p_val * q_val + r_val) / 2,
+            atol = 1e-5,
+        )
+        # Forward dx/dp
+        DiffOpt.set_forward_parameter(model, p, 1.0)
+        DiffOpt.set_forward_parameter(model, q, 0.0)
+        DiffOpt.set_forward_parameter(model, r, 0.0)
+        DiffOpt.forward_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_forward_variable(model, x),
+            -q_val / 2,
+            atol = 1e-5,
+        )
+        # Forward dx/dq
+        DiffOpt.set_forward_parameter(model, p, 0.0)
+        DiffOpt.set_forward_parameter(model, q, 1.0)
+        DiffOpt.forward_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_forward_variable(model, x),
+            -p_val / 2,
+            atol = 1e-5,
+        )
+        # Forward dx/dr
+        DiffOpt.set_forward_parameter(model, q, 0.0)
+        DiffOpt.set_forward_parameter(model, r, 1.0)
+        DiffOpt.forward_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_forward_variable(model, x),
+            -1 / 2,
+            atol = 1e-5,
+        )
+        # Forward combined direction: Δp=1, Δq=1, Δr=1
+        DiffOpt.set_forward_parameter(model, p, 1.0)
+        DiffOpt.set_forward_parameter(model, q, 1.0)
+        DiffOpt.set_forward_parameter(model, r, 1.0)
+        DiffOpt.forward_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_forward_variable(model, x),
+            -(q_val + p_val + 1) / 2,
+            atol = 1e-5,
+        )
+        # Reverse
+        DiffOpt.set_reverse_variable(model, x, 1.0)
+        DiffOpt.reverse_differentiate!(model)
+        @test isapprox(
+            DiffOpt.get_reverse_parameter(model, p),
+            -q_val / 2,
+            atol = 1e-5,
+        )
+        @test isapprox(
+            DiffOpt.get_reverse_parameter(model, q),
+            -p_val / 2,
+            atol = 1e-5,
+        )
+        @test isapprox(
+            DiffOpt.get_reverse_parameter(model, r),
+            -1 / 2,
+            atol = 1e-5,
+        )
+    end
+    return
+end
+
 end # module
 
 TestParameters.runtests()