diff --git a/src/chebyshev.jl b/src/chebyshev.jl
index 9c0f58b3..72bf4b72 100644
--- a/src/chebyshev.jl
+++ b/src/chebyshev.jl
@@ -12,6 +12,7 @@ mutable struct ChebyshevIterable{precT, matT, solT, vecT, realT <: Real}
     c::vecT
 
     α::realT
+    ζ::realT
 
     λ_avg::realT
     λ_diff::realT
@@ -32,17 +33,24 @@ function iterate(cheb::ChebyshevIterable, iteration::Int=start(cheb))
         return nothing
     end
 
-    T = eltype(cheb.x)
-
+    # Coupled two-term form of the Chebyshev recurrence. With μ = λ_avg, ρ = λ_diff
+    # and γ_k = T_k(μ / ρ), the residuals satisfy r_k = T_k((μI - A) / ρ) r_0 / γ_k.
+    # ζ holds the ratio γ_k / γ_{k-1}, which follows the Chebyshev three-term recurrence
+    # ζ_{k+1} = 2μ/ρ - 1/ζ_k, and the direction u_k = c_k + β_k u_{k-1} with β_k = 1/ζ_k².
     ldiv!(cheb.c, cheb.Pl, cheb.r)
 
-    if iteration == 1
-        cheb.α = T(2) / cheb.λ_avg
+    if iteration == 0
+        # First step is a Richardson step: r_1 = (I - A/μ) r_0.
+        cheb.ζ = cheb.λ_avg / cheb.λ_diff
+        cheb.α = inv(cheb.λ_avg)
         copyto!(cheb.u, cheb.c)
     else
-        β = (cheb.λ_diff * cheb.α / 2) ^ 2
-        cheb.α = inv(cheb.λ_avg - β)
+        # The three-term recurrence T_{k+1} = 2xT_k - T_{k-1} starts from T_{-1} = T_1,
+        # not T_0, so the second step uses half of the generic β.
+        β = iteration == 1 ? inv(2 * cheb.ζ^2) : inv(cheb.ζ^2)
         cheb.u .= cheb.c .+ β .* cheb.u
+        cheb.ζ = 2 * cheb.λ_avg / cheb.λ_diff - inv(cheb.ζ)
+        cheb.α = 2 / (cheb.λ_diff * cheb.ζ)
     end
 
     mul!(cheb.c, cheb.A, cheb.u)
@@ -84,7 +92,7 @@ function chebyshev_iterable!(x, A, b, λmin::Real, λmax::Real;
     tolerance = max(reltol * resnorm, abstol)
 
     ChebyshevIterable(Pl, A, x, r, u, c,
-        zero(real(T)),
+        zero(real(T)), zero(real(T)),
         λ_avg, λ_diff,
         resnorm, tolerance, maxiter, mv_products
     )
@@ -122,7 +130,7 @@ Solve Ax = b for symmetric, definite matrices A using Chebyshev iteration.
   tolerance for the stopping condition
   `|r_k| ≤ max(reltol * |r_0|, abstol)`, where `r_k = A * x_k - b`
   is the residual in the `k`th iteration;
-- `maxiter::Int = size(A, 2)`: maximum number of inner iterations of GMRES;
+- `maxiter::Int = size(A, 2)`: maximum number of iterations;
 - `Pl = Identity()`: left preconditioner;
 - `log::Bool = false`: keep track of the residual norm in each iteration;
 - `verbose::Bool = false`: print convergence information during the iterations.
diff --git a/test/chebyshev.jl b/test/chebyshev.jl
index 375ca00b..604984a4 100644
--- a/test/chebyshev.jl
+++ b/test/chebyshev.jl
@@ -66,6 +66,49 @@ Random.seed!(1234321)
     end
 end
 
+@testset "Residuals match the Chebyshev polynomial" begin
+    # r_k = T_k((μI - A) / ρ) r_0 / T_k(μ / ρ) with μ = (λ_max + λ_min) / 2 and
+    # ρ = (λ_max - λ_min) / 2. For a diagonal A (and diagonal Pl) this can be evaluated
+    # exactly per eigenvalue, so the residual history is checked against it directly.
+    function exact_resnorms(λ, r0, λ_min, λ_max, K)
+        μ = (λ_max + λ_min) / 2
+        ρ = (λ_max - λ_min) / 2
+        t = (μ .- λ) ./ ρ
+        s = μ / ρ
+        Tprev, Tcur = one.(t), copy(t)
+        Sprev, Scur = one(s), s
+        out = [norm(r0 .* Tcur ./ Scur)]
+        for k in 2:K
+            Tprev, Tcur = Tcur, 2 .* t .* Tcur .- Tprev
+            Sprev, Scur = Scur, 2 * s * Scur - Sprev
+            push!(out, norm(r0 .* Tcur ./ Scur))
+        end
+        out
+    end
+
+    n = 100
+    K = 40
+    for T in (Float32, Float64, ComplexF32, ComplexF64)
+        λ = collect(range(one(real(T)), real(T)(100), length = n))
+        A = Diagonal(T.(λ))
+        b = ones(T, n)
+        λ_min, λ_max = one(real(T)), real(T)(100)
+        rtol = 200 * eps(real(T))
+
+        x, history = chebyshev(A, b, λ_min, λ_max, maxiter = K, reltol = zero(real(T)), log = true)
+        @test history[:resnorm] ≈ exact_resnorms(λ, b, λ_min, λ_max, K) rtol = rtol
+
+        # Left preconditioner: the residual polynomial is applied to Pl \ A
+        p = collect(range(real(T)(0.5), real(T)(2), length = n))
+        Pl = Diagonal(T.(p))
+        λp = λ ./ p
+        λ_min, λ_max = minimum(λp), maximum(λp)
+        x, history = chebyshev(A, b, λ_min, λ_max, Pl = Pl, maxiter = K, reltol = zero(real(T)), log = true)
+        # `resnorm` is the unpreconditioned residual, which for diagonal Pl is p .* (Pl \ r)
+        @test history[:resnorm] ≈ exact_resnorms(λp, one.(λp), λ_min, λ_max, K) rtol = rtol
+    end
+end
+
 @testset "Termination criterion" begin
     for T in (Float32, Float64, ComplexF32, ComplexF64)
         A = T[ 2 -1  0