diff --git a/src/ADNLPProblems/chebyquad.jl b/src/ADNLPProblems/chebyquad.jl
new file mode 100644
index 00000000..e419942f
--- /dev/null
+++ b/src/ADNLPProblems/chebyquad.jl
@@ -0,0 +1,58 @@
+export chebyquad
+
+function chebyquad(; use_nls::Bool = false, kwargs...)
+  model = use_nls ? :nls : :nlp
+  return chebyquad(Val(model); kwargs...)
+end
+
+function Cheby(xj::Ti, i) where {Ti}
+  return (xj ≥ 1) * cosh(i * acosh(max(abs(xj), 1))) + (xj ≤ -1) * (-one(Ti))^(i) * cosh(i * acosh(max(abs(xj), 1))) + (abs(xj) ≤ 1) * cos(i * acos(min(max(xj, -1), 1)))
+end
+
+function chebyquad(
+  ::Val{:nlp};
+  n::Int = default_nvar,
+  m::Int = n,
+  type::Val{T} = Val(Float64),
+  chebyshev = Cheby,
+  kwargs...,
+) where {T}
+  m = max(m, n)
+  function f(x; n = length(x), m = m, chebyshev = chebyshev)
+    return 1 // 2 * sum((1 // n * sum( chebyshev(x[j], 2i) for j=1:n) + 1 // ((2i)^2 - 1))^2 for i = 1:Int(round(m / 2))) + 1 // 2 * sum((1 // n * sum( chebyshev(x[j], 2i - 1)  for j=1:n))^2 for i = 1:(Int(round(m / 2)) + mod(n, 2)))
+  end
+  x0 = T[j // (n + 1) for j=1:n]
+  return ADNLPModels.ADNLPModel(f, x0, name = "chebyquad"; kwargs...)
+end
+
+function chebyquad(
+  ::Val{:nls};
+  n::Int = default_nvar,
+  m::Int = n,
+  type::Val{T} = Val(Float64),
+  chebyshev = Cheby,
+  kwargs...,
+) where {T}
+  m = max(m, n)
+  function F!(r, x; n = length(x), m = length(r), chebyshev = chebyshev)
+    for i = 1:Int(round(m / 2))
+      r[2i] = 1 // n * sum(
+        (
+          chebyshev(x[j], 2i)
+        ) for j=1:n) + 1 // ((2i)^2 - 1)
+      r[2i - 1] = 1 // n * sum(
+        (
+          chebyshev(x[j], 2i - 1)
+        ) for j=1:n)
+    end
+    if mod(m, 2) == 1
+      r[m] = 1 // n * sum(
+        (
+          chebyshev(x[j], m)
+        ) for j=1:n)
+    end
+    return r
+  end
+  x0 = T[j // (n + 1) for j=1:n]
+  return ADNLPModels.ADNLSModel!(F!, x0, m, name = "chebyquad-nls"; kwargs...)
+end
diff --git a/src/Meta/chebyquad.jl b/src/Meta/chebyquad.jl
new file mode 100644
index 00000000..6a012f62
--- /dev/null
+++ b/src/Meta/chebyquad.jl
@@ -0,0 +1,26 @@
+chebyquad_meta = Dict(
+  :nvar => 100,
+  :variable_nvar => true,
+  :ncon => 0,
+  :variable_ncon => false,
+  :minimize => true,
+  :name => "chebyquad",
+  :has_equalities_only => false,
+  :has_inequalities_only => false,
+  :has_bounds => false,
+  :has_fixed_variables => false,
+  :objtype => :least_squares,
+  :contype => :unconstrained,
+  :best_known_lower_bound => -Inf,
+  :best_known_upper_bound => 500.0,
+  :is_feasible => true,
+  :defined_everywhere => missing,
+  :origin => :unknown,
+)
+get_chebyquad_nvar(; n::Integer = default_nvar, kwargs...) = n
+get_chebyquad_ncon(; n::Integer = default_nvar, kwargs...) = 0
+get_chebyquad_nlin(; n::Integer = default_nvar, kwargs...) = 0
+get_chebyquad_nnln(; n::Integer = default_nvar, kwargs...) = 0
+get_chebyquad_nequ(; n::Integer = default_nvar, kwargs...) = 0
+get_chebyquad_nineq(; n::Integer = default_nvar, kwargs...) = 0
+get_chebyquad_nls_nequ(; n::Integer = default_nvar, m::Int = n, kwargs...) = m
diff --git a/src/PureJuMP/chebyquad.jl b/src/PureJuMP/chebyquad.jl
new file mode 100644
index 00000000..faccb606
--- /dev/null
+++ b/src/PureJuMP/chebyquad.jl
@@ -0,0 +1,46 @@
+#
+#   The Chebychev quadrature problem in variable dimension, using the
+#   exact formula for the shifted Chebyshev polynomials.  This is a
+#   nonlinear least-squares problem with n groups. The Hessian is full.
+#
+#   Source: Problem 35 in
+#      J.J. More', B.S. Garbow and K.E. Hillstrom,
+#      "Testing Unconstrained Optimization Software",
+#      ACM Transactions on Mathematical Software, vol. 7(1), pp. 17-41, 1981.
+#   Also problem 58 in 
+#      A.R. Buckley,
+#      "Test functions for unconstrained minimization",
+#      TR 1989CS-3, Mathematics, statistics and computing centre,
+#      Dalhousie University, Halifax (CDN), 1989.
+#
+#   classification SBR2-AN-V-0
+
+export chebyquad
+
+function chebyquad(args...; n::Int = default_nvar, m::Int = n, kwargs...)
+  m = max(m, n)
+
+  nlp = Model()
+
+  x0 = [j/(n + 1) for j=1:n]
+  @variable(nlp, x[j = 1:n], start = x0[j])
+
+  # Chebyshev polynomial of the first kind, using explicit expression
+  @NLobjective(
+    nlp,
+    Min,
+    0.5 * sum((
+      1 / n * sum(
+        (
+          ifelse(x[j] ≥ 1, cosh(2i * acosh(x[j])), ifelse(x[j] ≤ -1, (-1)^(2i) * cosh(2i * acosh(-x[j])), cos(2i * acos(x[j]))))
+        ) for j=1:n) + 1 / ((2i)^2 - 1)
+    )^2 for i = 1:Int(round(m / 2))) + 0.5 * sum((
+      1 / n * sum(
+        (
+          ifelse(x[j] ≥ 1, cosh((2i - 1) * acosh(x[j])), ifelse(x[j] ≤ -1, (-1)^(2i - 1) * cosh((2i - 1) * acosh(-x[j])), cos((2i - 1) * acos(x[j]))))
+        ) for j=1:n)
+    )^2 for i = 1:(Int(round(m / 2)) + mod(n, 2)))
+  )
+
+  return nlp
+end