diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl
index 63e6745..682f5ad 100644
--- a/src/moi_nlp_model.jl
+++ b/src/moi_nlp_model.jl
@@ -4,6 +4,7 @@ mutable struct MathOptNLPModel <: AbstractNLPModel{Float64, Vector{Float64}}
   meta::NLPModelMeta{Float64, Vector{Float64}}
   eval::Union{MOI.AbstractNLPEvaluator, Nothing}
   lincon::LinearConstraints
+  quadcon::QuadraticConstraints
   obj::Objective
   counters::Counters
 end
@@ -33,7 +34,7 @@ function MathOptNLPModel(jmodel::JuMP.Model; hessian::Bool = true, name::String
     (nnln == 0 ? 0 : sum(length(nl_con.hess_I) for nl_con in eval.constraints)) : 0
 
   moimodel = backend(jmodel)
-  nlin, lincon, lin_lcon, lin_ucon = parser_MOI(moimodel)
+  nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon = parser_MOI(moimodel, nvar)
 
   if (eval ≠ nothing) && eval.has_nlobj
     obj = Objective("NONLINEAR", 0.0, spzeros(Float64, nvar), COO(), 0)
@@ -41,11 +42,11 @@ function MathOptNLPModel(jmodel::JuMP.Model; hessian::Bool = true, name::String
     obj = parser_objective_MOI(moimodel, nvar)
   end
 
-  ncon = nlin + nnln
-  lcon = vcat(lin_lcon, nl_lcon)
-  ucon = vcat(lin_ucon, nl_ucon)
-  nnzj = lincon.nnzj + nl_nnzj
-  nnzh = obj.nnzh + nl_nnzh
+  ncon = nlin + quadcon.nquad + nnln
+  lcon = vcat(lin_lcon, quad_lcon, nl_lcon)
+  ucon = vcat(lin_ucon, quad_ucon, nl_ucon)
+  nnzj = lincon.nnzj + quadcon.nnzj + nl_nnzj
+  nnzh = obj.nnzh + quadcon.nnzh + nl_nnzh
 
   meta = NLPModelMeta(
     nvar,
@@ -60,13 +61,13 @@ function MathOptNLPModel(jmodel::JuMP.Model; hessian::Bool = true, name::String
     nnzh = nnzh,
     lin = collect(1:nlin),
     lin_nnzj = lincon.nnzj,
-    nln_nnzj = nl_nnzj,
+    nln_nnzj = quadcon.nnzj + nl_nnzj,
     minimize = objective_sense(jmodel) == MOI.MIN_SENSE,
-    islp = (obj.type == "LINEAR") && (nnln == 0),
+    islp = (obj.type == "LINEAR") && (nnln == 0) && (quadcon.nquad == 0),
     name = name,
   )
 
-  return MathOptNLPModel(meta, eval, lincon, obj, Counters())
+  return MathOptNLPModel(meta, eval, lincon, quadcon, obj, Counters())
 end
 
 function NLPModels.obj(nlp::MathOptNLPModel, x::AbstractVector)
@@ -115,7 +116,13 @@ end
 
 function NLPModels.cons_nln!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVector)
   increment!(nlp, :neval_cons_nln)
-  MOI.eval_constraint(nlp.eval, c, x)
+  for i = 1:(nlp.quadcon.nquad)
+    qcon = nlp.quadcon[i]
+    c[i] = 0.5 * coo_sym_dot(qcon.hessian.rows, qcon.hessian.cols, qcon.hessian.vals, x, x) + dot(qcon.b, x)
+  end
+  if nlp.meta.nnln > nlp.quadcon.nquad
+    MOI.eval_constraint(nlp.eval, view(c, (nlp.quadcon.nquad + 1):(nlp.meta.nnln)), x)
+  end
   return c
 end
 
@@ -134,11 +141,16 @@ function NLPModels.jac_nln_structure!(
   rows::AbstractVector{<:Integer},
   cols::AbstractVector{<:Integer},
 )
-  jac_struct = MOI.jacobian_structure(nlp.eval)
-  for index = 1:(nlp.meta.nln_nnzj)
-    row, col = jac_struct[index]
-    rows[index] = row
-    cols[index] = col
+  quad_nnzj, jrows, jcols = nlp.quadcon.nnzj, nlp.quadcon.jrows, nlp.quadcon.jcols
+  rows[1:quad_nnzj] .= jrows
+  cols[1:quad_nnzj] .= jcols
+  if nlp.meta.nnln > nlp.quadcon.nquad
+    jac_struct = MOI.jacobian_structure(nlp.eval)
+    for index =  1:(nlp.meta.nln_nnzj - quad_nnzj)
+      row, col = jac_struct[index]
+      rows[quad_nnzj + index] = row + nlp.quadcon.nquad
+      cols[quad_nnzj + index] = col
+    end
   end
   return rows, cols
 end
@@ -151,7 +163,22 @@ end
 
 function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::AbstractVector)
   increment!(nlp, :neval_jac_nln)
-  MOI.eval_constraint_jacobian(nlp.eval, vals, x)
+  vals .= 0.0
+  k = 0
+  for i = 1:(nlp.quadcon.nquad)
+    qcon = nlp.quadcon[i]
+    for j=1:length(qcon.vec)
+      vals[k + j] = qcon.b[qcon.vec[j]]
+    end
+    nnzj = length(qcon.hessian.vals)
+    for i=1:nnzj
+      vals[k + i] += qcon.hessian.vals[i] * x[qcon.hessian.cols[i]]
+    end
+    k += nnzj
+  end
+  if nlp.meta.nnln > nlp.quadcon.nquad
+    MOI.eval_constraint_jacobian(nlp.eval, view(vals, (nlp.quadcon.nnzj + 1):(nlp.meta.nln_nnzj)), x)
+  end
   return vals
 end
 
@@ -305,10 +332,16 @@ function NLPModels.hess_structure!(
       cols[index] = nlp.obj.hessian.cols[index]
     end
   end
-  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0)
+  if nlp.quadcon.nquad > 0
+    for (index, tuple) in enumerate(nlp.quadcon.set)
+      rows[nlp.obj.nnzh + index] = tuple[1]
+      cols[nlp.obj.nnzh + index] = tuple[2]
+    end
+  end
+  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad)
     hesslag_struct = MOI.hessian_lagrangian_structure(nlp.eval)
-    for index = (nlp.obj.nnzh + 1):(nlp.meta.nnzh)
-      shift_index = index - nlp.obj.nnzh
+    for index = (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)
+      shift_index = index - nlp.obj.nnzh - nlp.quadcon.nnzh
       rows[index] = hesslag_struct[shift_index][1]
       cols[index] = hesslag_struct[shift_index][2]
     end
@@ -327,13 +360,23 @@ function NLPModels.hess_coord!(
   if nlp.obj.type == "QUADRATIC"
     vals[1:(nlp.obj.nnzh)] .= obj_weight .* nlp.obj.hessian.vals
   end
-  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0)
+  if nlp.quadcon.nquad > 0
+    quad_nnzh = nlp.quadcon.nnzh
+    k = 0
+    for i = 1:(nlp.quadcon.nquad)
+      qcon = nlp.quadcon[i]
+      nnzh = length(qcon.hessian.vals)
+      vals[(k + 1):(k + nnzh)] .= qcon.hessian.vals .* y[nlp.meta.nlin + i]
+      k += nnzh
+    end
+  end
+  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad)
     MOI.eval_hessian_lagrangian(
       nlp.eval,
-      view(vals, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)),
+      view(vals, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)),
       x,
       obj_weight,
-      view(y, nlp.meta.nln),
+      view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon))
     )
   end
   return vals
@@ -354,6 +397,7 @@ function NLPModels.hess_coord!(
     vals[(nlp.obj.nnzh + 1):(nlp.meta.nnzh)] .= 0.0
   end
   if nlp.obj.type == "NONLINEAR"
+    vals .= 0.0
     MOI.eval_hessian_lagrangian(nlp.eval, vals, x, obj_weight, zeros(nlp.meta.nnln))
   end
 
@@ -372,8 +416,18 @@ function NLPModels.hprod!(
   if (nlp.obj.type == "LINEAR") && (nlp.meta.nnln == 0)
     hv .= 0.0
   end
-  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0)
-    MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, view(y, nlp.meta.nln))
+  if nlp.quadcon.nquad > 0
+    for i = 1:(nlp.quadcon.nquad)
+      qcon = nlp.quadcon[i]
+      for k = 1:length(qcon.hessian.vals)
+        hv[qcon.hessian.rows[k]] += qcon.hessian.vals[k] * v[qcon.hessian.cols[k]]
+      end
+      hv[i] *= y[nlp.meta.nlin + i]
+    end
+  end
+  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad)
+    ind_nln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon)
+    MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, view(y, ind_nln))
   end
   if nlp.obj.type == "QUADRATIC"
     nlp.meta.nnln == 0 && (hv .= 0.0)
@@ -404,7 +458,8 @@ function NLPModels.hprod!(
     hv .*= obj_weight
   end
   if nlp.obj.type == "NONLINEAR"
-    MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, zeros(nlp.meta.nnln))
+    nnln = nlp.meta.nnln - nlp.quadcon.nquad
+    MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, zeros(nnln))
   end
   return hv
 end
diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl
index 29bf10b..01f1160 100644
--- a/src/moi_nls_model.jl
+++ b/src/moi_nls_model.jl
@@ -40,7 +40,7 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri
     (nnln == 0 ? 0 : sum(length(con.hess_I) for con in ceval.constraints)) : 0
 
   moimodel = backend(cmodel)
-  nlin, lincon, lin_lcon, lin_ucon = parser_MOI(moimodel)
+  nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon = parser_MOI(moimodel, nvar)
 
   nequ = nlinequ + nnlnequ
   Fnnzj = linequ.nnzj + nl_Fnnzj
diff --git a/src/utils.jl b/src/utils.jl
index 5910610..12913be 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -6,10 +6,18 @@ import NLPModels.increment!, NLPModels.decrement!
 using JuMP, MathOptInterface
 const MOI = MathOptInterface
 
+# VariableIndex
+const VI = MOI.VariableIndex  # VariableIndex(value)
+
 # ScalarAffineFunctions and VectorAffineFunctions
-const SAF = MOI.ScalarAffineFunction{Float64}
-const VAF = MOI.VectorAffineFunction{Float64}
-const AF = Union{SAF, VAF}
+const SAF = MOI.ScalarAffineFunction{Float64}  # ScalarAffineFunction{T}(terms, constant)
+const VAF = MOI.VectorAffineFunction{Float64}  # VectorAffineFunction{T}(terms, constants)
+const AF  = Union{SAF, VAF}
+
+# ScalarQuadraticFunctions and VectorQuadraticFunctions
+const SQF = MOI.ScalarQuadraticFunction{Float64}  # ScalarQuadraticFunction{T}(affine_terms, quadratic_terms, constant)
+const VQF = MOI.VectorQuadraticFunction{Float64}  # VectorQuadraticFunction{T}(affine_terms, quadratic_terms, constants)
+const QF  = Union{SQF, VQF}
 
 # AffLinSets and VecLinSets
 const ALS = Union{
@@ -21,10 +29,11 @@ const ALS = Union{
 const VLS = Union{MOI.Nonnegatives, MOI.Nonpositives, MOI.Zeros}
 const LS = Union{ALS, VLS}
 
+# Objective
 const VI = MOI.VariableIndex
-const SQF = MOI.ScalarQuadraticFunction{Float64}
 const OBJ = Union{VI, SAF, SQF}
 
+# Coordinate Matrix
 mutable struct COO
   rows::Vector{Int}
   cols::Vector{Int}
@@ -38,6 +47,25 @@ mutable struct LinearConstraints
   nnzj::Int
 end
 
+mutable struct QuadraticConstraint
+  hessian::COO
+  vec::Vector{Int}
+  b::SparseVector{Float64}
+end
+
+mutable struct QuadraticConstraints
+  qcons::Vector{QuadraticConstraint}
+  nquad::Int
+  nnzj::Int
+  jrows::Vector{Int}
+  jcols::Vector{Int}
+  nnzh::Int
+  set::Set{Tuple{Int,Int}}
+end
+
+Base.getindex(qcon::QuadraticConstraints, i::Integer) = qcon.qcons[i]
+Base.length(qcon::QuadraticConstraints) = qcon.nquad
+
 mutable struct LinearEquations
   jacobian::COO
   constants::Vector{Float64}
@@ -96,6 +124,48 @@ function coo_sym_dot(
   return xᵀAy
 end
 
+"""
+    jacobian_quad(qcons)
+
+`qcons` is a vector of `QuadraticConstraint` where each constraint has the form ½xᵀQᵢx + xᵀbᵢ.
+Compute the sparsity pattern of the jacobian [Q₁x + b₁; ...; Qₚx + bₚ]ᵀ of `qcons`.
+This function also allocates `qcons[i].vec`.
+"""
+function jacobian_quad(qcons)
+  jrows = Int[]
+  jcols = Int[]
+  nquad = length(qcons)
+  for i = 1 : nquad
+    # rows of Qᵢx + bᵢ with nonzeros coefficients
+    qcons[i].vec = unique(qcons[i].hessian.rows ∪ qcons[i].b.nzind)
+    for elt ∈ qcons[i].vec
+      push!(jcols, elt)
+      push!(jrows, i)
+    end
+  end
+  nnzj = length(jrows)
+  return nnzj, jrows, jcols
+end
+
+"""
+    hessian_quad(qcons)
+
+`qcons` is a vector of `QuadraticConstraint` where each constraint has the form ½xᵀQᵢx + xᵀbᵢ.
+Compute the sparsity pattern of the hessian ΣᵢQᵢ of `qcons`.
+"""
+function hessian_quad(qcons)
+  set = Set{Tuple{Int,Int}}()
+  nquad = length(qcons)
+  for i = 1 : nquad
+    con = qcons[i]
+    for tuple ∈ zip(con.hessian.rows, con.hessian.cols)
+      # Only disctinct tuples are stored in the set
+      push!(set, tuple)
+    end
+  end
+  return set
+end
+
 """
     parser_SAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon)
 
@@ -157,11 +227,64 @@ function parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_uco
 end
 
 """
-    parser_MOI(moimodel)
+    parser_SQF(fun, set, qcons, quad_lcon, quad_ucon)
 
-Parse linear constraints of a `MOI.ModelLike`.
+Parse a `ScalarQuadraticFunction` fun with its associated set.
+`qcons`, `quad_lcon`, `quad_ucon` are updated.
 """
-function parser_MOI(moimodel)
+function parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon)
+
+  b    = spzeros(Float64, nvar)
+  rows = Int[]
+  cols = Int[]
+  vals = Float64[]
+
+  # Parse a ScalarAffineTerm{Float64}(coefficient, variable_index)
+  for term in fun.affine_terms
+    b[term.variable.value] = term.coefficient
+  end
+
+  # Parse a ScalarQuadraticTerm{Float64}(coefficient, variable_index_1, variable_index_2)
+  for term in fun.quadratic_terms
+    i = term.variable_1.value
+    j = term.variable_2.value
+    if i ≥ j
+      push!(rows, i)
+      push!(cols, j)
+    else
+      push!(cols, j)
+      push!(rows, i)
+    end
+    push!(vals, term.coefficient)
+  end
+
+  if typeof(set) in (MOI.Interval{Float64}, MOI.GreaterThan{Float64})
+    push!(quad_lcon, -fun.constant + set.lower)
+  elseif typeof(set) == MOI.EqualTo{Float64}
+    push!(quad_lcon, -fun.constant + set.value)
+  else
+    push!(quad_lcon, -Inf)
+  end
+
+  if typeof(set) in (MOI.Interval{Float64}, MOI.LessThan{Float64})
+    push!(quad_ucon, -fun.constant + set.upper)
+  elseif typeof(set) == MOI.EqualTo{Float64}
+    push!(quad_ucon, -fun.constant + set.value)
+  else
+    push!(quad_ucon, Inf)
+  end
+
+  nnzh = length(vals)
+  qcon = QuadraticConstraint(COO(rows, cols, vals), Int[], b)
+  push!(qcons, qcon)
+end
+
+"""
+    parser_MOI(moimodel, nvar)
+
+Parse linear and quadratic constraints of a `MOI.ModelLike`.
+"""
+function parser_MOI(moimodel, nvar)
 
   # Variables associated to linear constraints
   nlin = 0
@@ -171,10 +294,16 @@ function parser_MOI(moimodel)
   lin_lcon = Float64[]
   lin_ucon = Float64[]
 
+  # Variables associated to quadratic constraints
+  nquad = 0
+  qcons = QuadraticConstraint[]
+  quad_lcon = Float64[]
+  quad_ucon = Float64[]
+
   contypes = MOI.get(moimodel, MOI.ListOfConstraintTypesPresent())
   for (F, S) in contypes
     F == VI && continue
-    F <: AF || @warn("Function $F is not supported.")
+    F <: AF || F <: SQF || @warn("Function $F is not supported.")
     S <: LS || @warn("Set $S is not supported.")
 
     conindices = MOI.get(moimodel, MOI.ListOfConstraintIndices{F, S}())
@@ -189,13 +318,22 @@ function parser_MOI(moimodel)
         parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon)
         nlin += set.dimension
       end
+      if typeof(fun) <: SQF
+        parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon)
+        nquad += 1
+      end
     end
   end
   coo = COO(linrows, lincols, linvals)
   nnzj = length(linvals)
   lincon = LinearConstraints(coo, nnzj)
 
-  return nlin, lincon, lin_lcon, lin_ucon
+  nnzj, jrows, jcols = jacobian_quad(qcons)
+  set = hessian_quad(qcons)
+  nnzh = length(set)
+  quadcon = QuadraticConstraints(qcons, nquad, nnzj, jrows, jcols, nnzh, set)
+
+  return nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon
 end
 
 """
diff --git a/test/nlp_problems/hs100.jl b/test/nlp_problems/hs100.jl
new file mode 100644
index 0000000..b4af448
--- /dev/null
+++ b/test/nlp_problems/hs100.jl
@@ -0,0 +1,24 @@
+function hs100(args...; kwargs...)
+  nlp = Model()
+  x0 = [1, 2, 0, 4, 0, 1, 1]
+  @variable(nlp, x[i = 1:7], start = x0[i])
+
+  @NLconstraint(nlp, 127 - 2 * x[1]^2 - 3 * x[2]^4 - x[3] - 4 * x[4]^2 - 5 * x[5] ≥ 0)
+  @constraint(nlp, 282 - 7 * x[1] - 3 * x[2] - 10 * x[3]^2 - x[4] + x[5] ≥ 0)
+  @constraint(nlp, 196 - 23 * x[1] - x[2]^2 - 6 * x[6]^2 + 8 * x[7] ≥ 0)
+  @constraint(nlp, -4 * x[1]^2 - x[2]^2 + 3 * x[1] * x[2] - 2 * x[3]^2 - 5 * x[6] + 11 * x[7] ≥ 0)
+
+  @NLobjective(
+    nlp,
+    Min,
+    (x[1] - 10)^2 +
+    5 * (x[2] - 12)^2 +
+    x[3]^4 +
+    3 * (x[4] - 11)^2 +
+    10 * x[5]^6 +
+    7 * x[6]^2 +
+    x[7]^4 - 4 * x[6] * x[7] - 10 * x[6] - 8 * x[7]
+  )
+
+  return nlp
+end
diff --git a/test/nlp_problems/hs61.jl b/test/nlp_problems/hs61.jl
new file mode 100644
index 0000000..5e422cc
--- /dev/null
+++ b/test/nlp_problems/hs61.jl
@@ -0,0 +1,12 @@
+"HS61 model"
+function hs61(args...; kwargs...)
+  nlp = Model()
+  @variable(nlp, x[i = 1:3], start = 0)
+
+  @constraint(nlp, 3 * x[1] - 2 * x[2]^2 - 7 == 0)
+  @constraint(nlp, 4 * x[1] - x[3]^2 - 11 == 0)
+
+  @NLobjective(nlp, Min, 4 * x[1]^2 + 2 * x[2]^2 + 2 * x[3]^2 - 33 * x[1] + 16 * x[2] - 24 * x[3])
+
+  return nlp
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 5306ee7..7a34589 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,7 +6,7 @@ nlp_problems = setdiff(NLPModelsTest.nlp_problems, ["MGH01Feas"])
 nls_problems = NLPModelsTest.nls_problems
 extra_nls_problems = ["HS30", "HS43", "MGH07", "nlsnohesspb"]
 
-for problem in lowercase.(nlp_problems ∪ ["nohesspb"])
+for problem in lowercase.(nlp_problems ∪ ["nohesspb", "hs61", "hs100"])
   include(joinpath("nlp_problems", "$problem.jl"))
 end
 
diff --git a/test/test_moi_nlp_model.jl b/test/test_moi_nlp_model.jl
index 162694d..cd4bb20 100644
--- a/test/test_moi_nlp_model.jl
+++ b/test/test_moi_nlp_model.jl
@@ -11,7 +11,7 @@ println("Testing MathOptNLPModel")
   "‖c(x₀)‖"
 )
 # Test that every problem can be instantiated.
-for prob in Symbol.(lowercase.(nlp_problems ∪ ["nohesspb"]))
+for prob in Symbol.(lowercase.(nlp_problems ∪ ["nohesspb", "hs61", "hs100"]))
   prob_fn = eval(prob)
   nlp = MathOptNLPModel(prob_fn(), hessian = (prob != :nohesspb))
   n = nlp.meta.nvar
@@ -23,3 +23,47 @@ for prob in Symbol.(lowercase.(nlp_problems ∪ ["nohesspb"]))
   @printf("%-15s  %4d  %4d  %10.4e  %10.4e  %10s\n", prob, n, m, fx, ngx, ncx)
 end
 println()
+
+function hs219(args...; kwargs...)
+  nlp = Model()
+  x0 = [10, 10, 10, 10]
+  @variable(nlp, x[i = 1:4], start = x0[i])
+
+  @constraint(nlp, x[1]^2 - x[2] - x[4]^2 == 0)
+  @NLconstraint(nlp, x[2] - x[1]^3 - x[3]^2 == 0)
+
+  @NLobjective(
+    nlp,
+    Min,
+    -x[1]
+  )
+
+  return nlp
+end
+
+@testset "Testing quadratic constraints with JuMP" begin
+  g(x) = [-1., 0., 0., 0.]
+  Hess(x) = zeros(4, 4)
+  function Hess(x, y) 
+    H = zeros(4, 4)
+    H[1, 1] = 2 * y[1] - y[2] * 6 * x[1]
+    H[3, 3] = - 2 * y[2]
+    H[4, 4] = - 2 * y[1]
+    return H
+  end
+  J(x) = [
+      2x[1] -1 0 -2x[4];
+      -3x[1]^2 1 -2x[3] 0
+  ]
+
+  jump = hs219()
+  nlp = MathOptNLPModel(jump)
+  x1 = nlp.meta.x0
+  y1 = ones(nlp.meta.ncon)
+  v1 = 2 * ones(nlp.meta.nvar)
+  @test jac(nlp, x1) ≈ J(x1)
+  @test hess(nlp, x1) ≈ Hess(x1)
+  @test hess(nlp, x1, y1) ≈ Hess(x1, y1)
+  @test hprod(nlp, x1, x1) ≈ Hess(x1) * x1
+  @test hprod(nlp, x1, y1, v1) ≈ Hess(x1, y1) * v1
+end