diff --git a/CHANGELOG.md b/CHANGELOG.md
index dd9c7964..f2fb8fc3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,10 @@
 
 ### Fixed
 
+- A `POLY` section's cl and cd no longer jump at |α| = 20°: they blend from the
+  α-polynomial into the flat-plate post-stall form over 20°–25°. Past 20°, cl now has
+  the sign of α and cd stays positive, and cm holds its ±20° value instead of running
+  its quadratic unbounded. Results inside ±20° are unchanged.
 - The docs define the body frame as KiteUtils' `KA` frame: x from LE to TE, y towards
   the right tip, z = x × y up. The separate `KB` definition, whose `Z = Y × X` pointed
   down, is gone; the solver's frame is unchanged. CL, CD and CS are documented as
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index f7482533..5e889e49 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -124,6 +124,9 @@ reads_from
 same_knots
 polar_knots
 window_alpha
+blend_post_stall
+flat_plate_cl
+flat_plate_cd
 assemble_polar_matrix
 load_polar_data
 load_matrix_polar_data
diff --git a/docs/src/private_types.md b/docs/src/private_types.md
index 0941c427..3b6d810a 100644
--- a/docs/src/private_types.md
+++ b/docs/src/private_types.md
@@ -19,6 +19,8 @@ SemiInfiniteFilament
 ```@docs
 LEI_AIRFOIL_BREUKELS
 SCAN_KNOT_MAX
+POLY_FIT_LIMIT
+POLY_BLEND_WIDTH
 ```
 
 ### Airfoil aerodynamics (AirfoilAero)
diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl
index 6f473f1c..8fcf6075 100644
--- a/src/VortexStepMethod.jl
+++ b/src/VortexStepMethod.jl
@@ -254,7 +254,9 @@ Enumeration of the implemented aerodynamic models. See also: [`AeroData`](@ref)
 
 # Elements
 - `POLY`: α-polynomial coefficients for cl/cd/cm (e.g. Breukels LEI coeffs, generated
-  by the `AirfoilAero` package). Core only evaluates the polynomial.
+  by the `AirfoilAero` package), fitted inside ±20°. Past that, cl and cd blend into the
+  flat-plate forms 2 sin α |sin α| cos α and 2 |sin α|³ by ±25°, and cm holds its ±20°
+  value.
 - `POLAR_VECTORS`: Polar vectors as function of alpha (lookup tables with interpolation).
   A panel's table may be rewritten at run time by [`set_polar!`](@ref
   VortexStepMethod.set_polar!), which is how a live polar source regenerates it from the
diff --git a/src/panel.jl b/src/panel.jl
index 3a211857..34420d23 100644
--- a/src/panel.jl
+++ b/src/panel.jl
@@ -489,6 +489,32 @@ function calculate_relative_alpha_and_velocity(panel::Panel, induced_velocity)
     return flow.alpha, flow.v_eff
 end
 
+"Edge of the `POLY` α-polynomial fit [deg]."
+const POLY_FIT_LIMIT = 20.0
+"Width over which a `POLY` cl and cd blend into the flat-plate post-stall form [deg]."
+const POLY_BLEND_WIDTH = 5.0
+
+"Flat-plate post-stall lift coefficient at `alpha` [rad]."
+flat_plate_cl(alpha) = 2 * cos(alpha) * sin(alpha) * abs(sin(alpha))
+"Flat-plate post-stall drag coefficient at `alpha` [rad]."
+flat_plate_cd(alpha) = 2 * abs(sin(alpha))^3
+
+"""
+    blend_post_stall(coeffs, flat_plate, alpha)
+
+`POLY` coefficient at `alpha` [rad]: the α-polynomial `coeffs` (α in degrees) inside
+±`POLY_FIT_LIMIT`, the flat-plate form `flat_plate(alpha)` past ±(`POLY_FIT_LIMIT` +
+`POLY_BLEND_WIDTH`), and a cubic smoothstep between the two.
+"""
+function blend_post_stall(coeffs, flat_plate, alpha)
+    alpha_deg = rad2deg(alpha)
+    polynomial = evalpoly(alpha_deg, coeffs)
+    abs(alpha_deg) <= POLY_FIT_LIMIT && return polynomial
+    progress = min((abs(alpha_deg) - POLY_FIT_LIMIT) / POLY_BLEND_WIDTH, 1)
+    weight = progress^2 * (3 - 2progress)
+    return (1 - weight) * polynomial + weight * flat_plate(alpha)
+end
+
 """
     calculate_cl(panel::Panel, alpha)
     calculate_cl(panel::Panel, alpha, delta)
@@ -507,11 +533,7 @@ function calculate_cl(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td}
     isnan(alpha) && return R(NaN)
     alpha = window_alpha(panel, alpha)
     if panel.aero_model == POLY
-        cl = evalpoly(rad2deg(alpha), panel.cl_coeffs)
-        if abs(alpha) > (π/9)
-            cl = 2 * cos(alpha) * sin(alpha)^2
-        end
-        return R(cl)
+        return R(blend_post_stall(panel.cl_coeffs, flat_plate_cl, alpha))
     elseif panel.aero_model == INVISCID
         return R(2π * alpha)
     end
@@ -537,10 +559,7 @@ function calculate_cd(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td}
     isnan(alpha) && return R(NaN)
     alpha = window_alpha(panel, alpha)
     if panel.aero_model == POLY
-        if abs(alpha) > (π/9)  # Outside ±20 degrees
-            return R(2 * sin(alpha)^3)
-        end
-        return R(evalpoly(rad2deg(alpha), panel.cd_coeffs))
+        return R(blend_post_stall(panel.cd_coeffs, flat_plate_cd, alpha))
     elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES)
         cd_interp = panel.cd_interp
         cd_interp === nothing &&
@@ -567,7 +586,8 @@ function calculate_cm(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td}
     isnan(alpha) && return R(NaN)
     alpha = window_alpha(panel, alpha)
     if panel.aero_model == POLY
-        return R(evalpoly(rad2deg(alpha), panel.cm_coeffs))
+        return R(evalpoly(clamp(rad2deg(alpha), -POLY_FIT_LIMIT, POLY_FIT_LIMIT),
+                          panel.cm_coeffs))
     elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES)
         cm_interp = panel.cm_interp
         cm_interp === nothing &&
diff --git a/test/panel/test_panel.jl b/test/panel/test_panel.jl
index 63fc8231..3f5c237d 100644
--- a/test/panel/test_panel.jl
+++ b/test/panel/test_panel.jl
@@ -1,4 +1,6 @@
 using VortexStepMethod: Panel, Section, calculate_relative_alpha_and_relative_velocity, calculate_cl, calculate_cd_cm, reinit!, INVISCID, POLAR_VECTORS, MVec3
+using VortexStepMethod: POLY, calculate_cd, calculate_cm
+using VortexStepMethod.AirfoilAero: lei_poly_coeffs
 using VortexStepMethod: panel_axes, panel_span_vector
 using Interpolations: linear_interpolation, Line
 using LinearAlgebra
@@ -216,3 +218,41 @@ end
         @test isapprox(flipped.x_airf, axes.x_airf; atol=1e-12)
     end
 end
+
+@testset "POLY polar is continuous past its ±20° fit range" begin
+    coeffs = lei_poly_coeffs(0.1, 0.08)
+    panel = create_panel(Section([0.0, 1.0, 0.0], [1.0, 1.0, 0.0], POLY, coeffs),
+                         Section([0.0, -1.0, 0.0], [1.0, -1.0, 0.0], POLY, coeffs))
+    polar(alpha) = (calculate_cl(panel, alpha), calculate_cd(panel, alpha),
+                    calculate_cm(panel, alpha))
+    alphas = deg2rad.(-60:0.01:60)
+
+    @testset "the polynomial is untouched inside ±20°" begin
+        for alpha in deg2rad.(-20:0.5:20)
+            @test polar(alpha) == evalpoly.(rad2deg(alpha), coeffs)
+        end
+    end
+    @testset "no coefficient jumps between samples 0.01° apart" begin
+        values = reduce(hcat, collect.(polar.(alphas)))
+        @test maximum(abs, diff(values; dims=2)) < 0.01
+    end
+    @testset "cl and cd slopes do not kink at the blend edges" begin
+        for coefficient in (calculate_cl, calculate_cd)
+            values = coefficient.(Ref(panel), alphas)
+            @test maximum(abs, diff(diff(values))) < 1e-4
+        end
+    end
+    @testset "post-stall lift is odd and drag even and positive" begin
+        for alpha in deg2rad.(25:5:85)
+            @test calculate_cl(panel, -alpha) ≈ -calculate_cl(panel, alpha)
+            @test calculate_cd(panel, -alpha) ≈ calculate_cd(panel, alpha)
+            @test calculate_cd(panel, alpha) > 0
+        end
+    end
+    @testset "cm holds its ±20° value past the fit range" begin
+        for alpha_deg in 20:5:90, side in (-1, 1)
+            @test calculate_cm(panel, deg2rad(side * alpha_deg)) ≈
+                  evalpoly(side * 20.0, coeffs[3])
+        end
+    end
+end