diff --git a/src/ADNLPProblems/chebyquad.jl b/src/ADNLPProblems/chebyquad.jl
new file mode 100644
index 00000000..a5e8e445
--- /dev/null
+++ b/src/ADNLPProblems/chebyquad.jl
@@ -0,0 +1,76 @@
+export chebyquad
+
+function chebyquad(; use_nls::Bool = false, kwargs...)
+  model = use_nls ? :nls : :nlp
+  return chebyquad(Val(model); kwargs...)
+end
+
+function Cheby(xj, i::Integer)
+  # Evaluate T_i(x) via recurrence to avoid domain-branching and AD/tracer issues.
+  if i == 0
+    return one(xj)
+  elseif i == 1
+    return xj
+  end
+
+  tk_minus_1 = one(xj)
+  tk = xj
+  for _ = 2:i
+    tk_plus_1 = 2 * xj * tk - tk_minus_1
+    tk_minus_1 = tk
+    tk = tk_plus_1
+  end
+  return tk
+end
+
+function chebyquad(
+  ::Val{:nlp};
+  n::Int = default_nvar,
+  m::Int = n,
+  type::Type{T} = Float64,
+  chebyshev = Cheby,
+  kwargs...,
+) where {T}
+  m = max(m, n)
+  function f(x; n = length(x), m = m, chebyshev = chebyshev)
+    return (one(T) / 2) * sum(
+      ((one(T) / n) * sum(chebyshev(x[j], 2i) for j = 1:n) + one(T) / ((2i)^2 - 1))^2 for
+      i = 1:div(m, 2)
+    ) +
+           (one(T) / 2) * sum(
+      ((one(T) / n) * sum(chebyshev(x[j], 2i - 1) for j = 1:n))^2 for i = 1:div(m + 1, 2)
+    )
+  end
+  x0 = Vector{T}(undef, n)
+  for j = 1:n
+    x0[j] = j * one(T) / (n + one(T))
+  end
+  return ADNLPModels.ADNLPModel(f, x0, name = "chebyquad"; kwargs...)
+end
+
+function chebyquad(
+  ::Val{:nls};
+  n::Int = default_nvar,
+  m::Int = n,
+  type::Type{T} = 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:div(m, 2)
+      r[2i] = (one(T) / n) * sum(chebyshev(x[j], 2i) for j = 1:n) + one(T) / ((2i)^2 - 1)
+      r[2i - 1] = (one(T) / n) * sum(chebyshev(x[j], 2i - 1) for j = 1:n)
+    end
+    if mod(m, 2) == 1
+      r[m] = (one(T) / n) * sum(chebyshev(x[j], m) for j = 1:n)
+    end
+    return r
+  end
+  x0 = Vector{T}(undef, n)
+  for j = 1:n
+    x0[j] = j * one(T) / (n + one(T))
+  end
+  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..65b36e39
--- /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...) = max(m, n)
diff --git a/src/PureJuMP/chebyquad.jl b/src/PureJuMP/chebyquad.jl
new file mode 100644
index 00000000..248a4c80
--- /dev/null
+++ b/src/PureJuMP/chebyquad.jl
@@ -0,0 +1,59 @@
+#
+#   The Chebyshev 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(
+            2 * x[j] - 1 ≥ 1,
+            cosh(2i * acosh(2 * x[j] - 1)),
+            ifelse(
+              2 * x[j] - 1 ≤ -1,
+              (-1)^(2i) * cosh(2i * acosh(1 - 2 * x[j])),
+              cos(2i * acos(2 * x[j] - 1)),
+            ),
+          ) for j = 1:n
+        ) + 1 / ((2i)^2 - 1)
+      )^2 for i = 1:div(m, 2)
+    ) +
+    0.5 * sum(
+      (
+        1 / n * sum(
+          ifelse(
+            2 * x[j] - 1 ≥ 1,
+            cosh((2i - 1) * acosh(2 * x[j] - 1)),
+            ifelse(
+              2 * x[j] - 1 ≤ -1,
+              (-1)^(2i - 1) * cosh((2i - 1) * acosh(1 - 2 * x[j])),
+              cos((2i - 1) * acos(2 * x[j] - 1)),
+            ),
+          ) for j = 1:n
+        )
+      )^2 for i = 1:div(m + 1, 2)
+    )
+  )
+  return nlp
+end