diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4ca2653a..e08187bc 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,9 @@
 
 ### Added
 
+- `stability_derivatives` gives the force and moment coefficients and their derivatives
+  with respect to angle of attack and sideslip, and `trim_angle` the angles of attack at
+  which `CMy` changes sign, with the slope that says whether each trim is stable.
 - Spanwise-flow viscous drag correction (Gaunaa et al. 2024,
   doi:10.1088/1742-6596/2767/2/022068): each section gets a drag increment and a force
   along its span from the flow across it, in `solve!`, `solve` and `linearize`. Opt-in
diff --git a/docs/src/functions.md b/docs/src/functions.md
index fb805d29..76d8dc8a 100644
--- a/docs/src/functions.md
+++ b/docs/src/functions.md
@@ -102,6 +102,8 @@ solve!
 solve_base!
 reinit!(body_aero::BodyAerodynamics{P, W, T}) where {P, W, T}
 linearize
+stability_derivatives
+trim_angle
 calculate_results
 ```
 
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index 15f636ef..c508a6f4 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -25,6 +25,10 @@ calculate_cd
 calculate_cm
 calculate_cd_cm
 set_pitch_rate_dist!
+apparent_wind
+coeffs_at_angles
+nose_down
+bisect_sign_change
 calculate_relative_alpha_and_velocity
 calculate_relative_alpha_and_relative_velocity
 update_effective_angle_of_attack!
diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl
index 0521ddfc..ab30ae28 100644
--- a/src/VortexStepMethod.jl
+++ b/src/VortexStepMethod.jl
@@ -32,6 +32,7 @@ export slice_args, preview_args
 export ObjWing, Section, Wing, refine!, reinit!
 export BodyAerodynamics
 export Solver, VSMSolution, linearize, solve, solve!, solve_base!, calc_forces!
+export stability_derivatives, trim_angle
 export SolveFailure
 export calculate_results
 export add_section!, set_va!, section_pitch_rate
@@ -432,6 +433,7 @@ include("panel.jl")
 include("body_aerodynamics.jl")
 include("wake.jl")
 include("solver.jl")
+include("stability.jl")
 
 include("plotting_helpers.jl")
 
diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index edac27a7..8a64ce89 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -1135,24 +1135,19 @@ function set_va!(body_aero::BodyAerodynamics, va_vec_dist::AbstractMatrix;
 end
 
 """
-    set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
-
-Set velocity array from VSM settings configuration.
+    apparent_wind(alpha, beta, wind_speed)
 
-This convenience method extracts flight conditions from VSMSettings and 
-constructs the velocity vector in the body reference frame based on:
-- Wind speed from settings.condition.wind_speed
-- Angle of attack from settings.condition.alpha (converted from degrees)
-- Sideslip angle from settings.condition.beta (converted from degrees)
+Apparent wind vector in the body frame [m/s] at angle of attack `alpha` [rad], sideslip
+`beta` [rad] and `wind_speed` [m/s].
+"""
+apparent_wind(alpha, beta, wind_speed) =
+    wind_speed .* [cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)]
 
-The velocity vector is constructed as:
-- X_b (forward): wind_speed * cos(α) * cos(β)  
-- Y_b (right): wind_speed * sin(β)
-- Z_b (down): wind_speed * sin(α) * cos(β)
+"""
+    set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
 
-# Arguments
-- `body_aero::BodyAerodynamics`: The aerodynamic body to modify
-- `settings::VSMSettings`: Settings object containing flight conditions
+Set the uniform inflow of `body_aero` to the [`apparent_wind`](@ref) at the `alpha` and
+`beta` [°] and `wind_speed` [m/s] of `settings.condition`.
 
 # Example
 ```julia
@@ -1162,15 +1157,8 @@ set_va!(body_aero, settings)
 ```
 """
 function set_va!(body_aero::BodyAerodynamics, settings::VSMSettings)
-    α = deg2rad(settings.condition.alpha)
-    β = deg2rad(settings.condition.beta)
-    wind_speed = settings.condition.wind_speed
-    
-    va_vec = wind_speed * [
-        cos(α)*cos(β),  # X_b (forward)
-        sin(β),         # Y_b (right)
-        sin(α)*cos(β)   # Z_b (down)
-    ]
-    
+    condition = settings.condition
+    va_vec = apparent_wind(deg2rad(condition.alpha), deg2rad(condition.beta),
+        condition.wind_speed)
     set_va!(body_aero, va_vec)
 end
diff --git a/src/stability.jl b/src/stability.jl
new file mode 100644
index 00000000..9ca803cf
--- /dev/null
+++ b/src/stability.jl
@@ -0,0 +1,92 @@
+"""
+    stability_derivatives(solver, body_aero, alpha, beta, wind_speed; kwargs...)
+
+Aerodynamic coefficients `[CFx, CFy, CFz, CMx, CMy, CMz]` of `body_aero` at angle of attack
+`alpha` [rad], sideslip `beta` [rad] and `wind_speed` [m/s], and their derivatives with
+respect to `alpha` and `beta` [1/rad], at the rotation rate `body_aero.omega` and with
+moments about `solver.reference_point`. `kwargs` go to [`linearize`](@ref), which leaves
+`body_aero` at this inflow.
+
+Returns `(coeffs, dalpha, dbeta, converged)`.
+"""
+function stability_derivatives(solver::Solver, body_aero::BodyAerodynamics, alpha, beta,
+        wind_speed; kwargs...)
+    va_vec = apparent_wind(alpha, beta, wind_speed)
+    jac, results, converged = linearize(solver, body_aero, va_vec;
+        theta_idxs=nothing, va_idxs=1:3, aero_coeffs=true, kwargs...)
+    dva_dalpha = ForwardDiff.derivative(
+        angle -> apparent_wind(angle, beta, wind_speed), alpha)
+    dva_dbeta = ForwardDiff.derivative(
+        angle -> apparent_wind(alpha, angle, wind_speed), beta)
+    coeff_jac = jac[1:6, :]
+    return (coeffs=results[1:6], dalpha=coeff_jac * dva_dalpha,
+        dbeta=coeff_jac * dva_dbeta, converged)
+end
+
+"""
+    trim_angle(solver, body_aero, beta, wind_speed; alpha_range=deg2rad.(-5:2:15),
+               alpha_tol=1e-5, backend=AutoForwardDiff())
+
+Angles of attack [rad] at which `CMy` of `body_aero` about `solver.reference_point` changes
+sign between neighbouring entries of `alpha_range`, bisected to `alpha_tol` [rad], at
+sideslip `beta` [rad] and `wind_speed` [m/s]. Returns one `(alpha, dCMy_dalpha)` per trim,
+the slope [1/rad] from [`stability_derivatives`](@ref) with `backend`; a trim is statically
+stable where `dCMy_dalpha < 0`. Throws a [`SolveFailure`](@ref) if a solve misses the
+solver's tolerances.
+"""
+function trim_angle(solver::Solver, body_aero::BodyAerodynamics, beta, wind_speed;
+        alpha_range=deg2rad.(-5:2:15), alpha_tol=1e-5, backend=AutoForwardDiff())
+    is_nose_down = alpha -> nose_down(solver, body_aero, alpha, beta, wind_speed)
+    nose_down_range = is_nose_down.(alpha_range)
+    trims = @NamedTuple{alpha::Float64, dCMy_dalpha::Float64}[]
+    for i in 1:length(alpha_range)-1
+        nose_down_range[i] == nose_down_range[i+1] && continue
+        alpha = bisect_sign_change(is_nose_down, alpha_range[i], alpha_range[i+1],
+            alpha_tol)
+        derivatives = stability_derivatives(solver, body_aero, alpha, beta, wind_speed;
+            backend, throw_on_fail=true)
+        push!(trims, (alpha=alpha, dCMy_dalpha=derivatives.dalpha[5]))
+    end
+    return trims
+end
+
+"""
+    coeffs_at_angles(solver, body_aero, alpha, beta, wind_speed)
+
+Aerodynamic coefficients `[CFx, CFy, CFz, CMx, CMy, CMz]` of `body_aero` solved at angle of
+attack `alpha` [rad], sideslip `beta` [rad] and `wind_speed` [m/s], at the rotation rate
+`body_aero.omega`. Throws a [`SolveFailure`](@ref) if the solve misses the solver's
+tolerances.
+"""
+function coeffs_at_angles(solver, body_aero, alpha, beta, wind_speed)
+    set_va!(body_aero, apparent_wind(alpha, beta, wind_speed), body_aero.omega)
+    sol = solve!(solver, body_aero; throw_on_fail=true)
+    return [sol.force_coeffs; sol.moment_coeffs]
+end
+
+"""
+    nose_down(solver, body_aero, alpha, beta, wind_speed)
+
+Whether `CMy` from [`coeffs_at_angles`](@ref) is negative.
+"""
+nose_down(solver, body_aero, alpha, beta, wind_speed) =
+    coeffs_at_angles(solver, body_aero, alpha, beta, wind_speed)[5] < 0
+
+"""
+    bisect_sign_change(predicate, low, high, tol)
+
+Bisect `[low, high]`, across which the boolean `predicate` flips, to a width of `tol` and
+return the midpoint.
+"""
+function bisect_sign_change(predicate, low, high, tol)
+    predicate_low = predicate(low)
+    while high - low > tol
+        middle = (low + high) / 2
+        if predicate(middle) == predicate_low
+            low = middle
+        else
+            high = middle
+        end
+    end
+    return (low + high) / 2
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index b77976d4..3077e22b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -58,6 +58,7 @@ function include_selected_tests()
     should_run_test("solver/test_forwarddiff.jl") && include("solver/test_forwarddiff.jl")
     should_run_test("solver/test_backend_comparison.jl") && include("solver/test_backend_comparison.jl")
     should_run_test("solver/test_unrefined_dist.jl") && include("solver/test_unrefined_dist.jl")
+    should_run_test("solver/test_stability.jl") && include("solver/test_stability.jl")
     should_run_test("verification/test_verification.jl") && include("verification/test_verification.jl")
     should_run_test("VortexStepMethod/test_VortexStepMethod.jl") && include("VortexStepMethod/test_VortexStepMethod.jl")
     should_run_test("wake/test_wake.jl") && include("wake/test_wake.jl")
diff --git a/test/solver/test_stability.jl b/test/solver/test_stability.jl
new file mode 100644
index 00000000..3e26529f
--- /dev/null
+++ b/test/solver/test_stability.jl
@@ -0,0 +1,92 @@
+using VortexStepMethod
+using VortexStepMethod: coeffs_at_angles
+using Test
+
+"""
+    trimmable_wing_aero(cm)
+
+A rectangular wing with lift slope 2π, `cl = 0.1` at zero angle of attack and a constant
+section `cm`.
+"""
+function trimmable_wing_aero(cm)
+    chord, span = 1.0, 6.0
+    alpha_range = deg2rad.(-10.0:5.0:20.0)
+    polar = (alpha_range, 2π .* alpha_range .+ 0.1, fill(0.02, length(alpha_range)),
+             fill(cm, length(alpha_range)))
+    wing = Wing(10)
+    add_section!(wing, [0.0, span / 2, 0.0], [chord, span / 2, 0.0], POLAR_VECTORS, polar)
+    add_section!(wing, [0.0, -span / 2, 0.0], [chord, -span / 2, 0.0], POLAR_VECTORS,
+                 polar)
+    refine!(wing)
+    return BodyAerodynamics([wing])
+end
+
+@testset "stability_derivatives match central differences of solve!" begin
+    body_aero = trimmable_wing_aero(0.05)
+    solver = Solver(body_aero; reference_point=[0.25, 0.5, 0.1], use_gamma_prev=false)
+    alpha, beta, wind_speed, step = deg2rad(4.0), deg2rad(3.0), 20.0, 1e-4
+
+    derivatives = stability_derivatives(solver, body_aero, alpha, beta, wind_speed)
+    @test derivatives.converged
+    @test derivatives.coeffs ≈
+          coeffs_at_angles(solver, body_aero, alpha, beta, wind_speed)
+
+    central_difference(coeffs_plus, coeffs_minus) = (coeffs_plus - coeffs_minus) / 2step
+    dalpha = central_difference(
+        coeffs_at_angles(solver, body_aero, alpha + step, beta, wind_speed),
+        coeffs_at_angles(solver, body_aero, alpha - step, beta, wind_speed))
+    dbeta = central_difference(
+        coeffs_at_angles(solver, body_aero, alpha, beta + step, wind_speed),
+        coeffs_at_angles(solver, body_aero, alpha, beta - step, wind_speed))
+    @test !iszero(dbeta)
+    @test derivatives.dalpha ≈ dalpha rtol = 1e-4 atol = 1e-6
+    @test derivatives.dbeta ≈ dbeta rtol = 1e-4 atol = 1e-6
+end
+
+@testset "trim_angle finds where CMy changes sign" begin
+    beta, wind_speed = 0.0, 20.0
+
+    @testset "moments about the leading edge: stable trim" begin
+        body_aero = trimmable_wing_aero(0.05)
+        solver = Solver(body_aero)
+        trims = trim_angle(solver, body_aero, beta, wind_speed)
+        @test length(trims) == 1
+        trim = only(trims)
+        trim_coeffs = coeffs_at_angles(solver, body_aero, trim.alpha, beta, wind_speed)
+        @test abs(trim_coeffs[5]) < 1e-5
+        @test trim.dCMy_dalpha < 0
+        derivatives = stability_derivatives(solver, body_aero, trim.alpha, beta, wind_speed)
+        @test trim.dCMy_dalpha ≈ derivatives.dalpha[5]
+    end
+
+    @testset "moments about the trailing edge: unstable trim" begin
+        body_aero = trimmable_wing_aero(-0.05)
+        solver = Solver(body_aero; reference_point=[1.0, 0.0, 0.0])
+        trim = only(trim_angle(solver, body_aero, beta, wind_speed))
+        trim_coeffs = coeffs_at_angles(solver, body_aero, trim.alpha, beta, wind_speed)
+        @test abs(trim_coeffs[5]) < 1e-5
+        @test trim.dCMy_dalpha > 0
+    end
+
+    @testset "a NONLIN solver with backend=nothing finds the same trim" begin
+        body_aero = trimmable_wing_aero(0.05)
+        trim_loop = only(trim_angle(Solver(body_aero), body_aero, beta, wind_speed))
+        solver = Solver(body_aero; solver_type=NONLIN)
+        trim = only(trim_angle(solver, body_aero, beta, wind_speed; backend=nothing))
+        @test trim.alpha ≈ trim_loop.alpha atol = 1e-4
+        @test trim.dCMy_dalpha ≈ trim_loop.dCMy_dalpha rtol = 1e-4
+    end
+
+    @testset "a solve that misses the tolerances throws" begin
+        body_aero = trimmable_wing_aero(0.05)
+        solver = Solver(body_aero; max_iterations=2)
+        @test_throws SolveFailure trim_angle(solver, body_aero, beta, wind_speed)
+    end
+
+    @testset "no sign change in alpha_range: no trim" begin
+        body_aero = trimmable_wing_aero(0.05)
+        solver = Solver(body_aero)
+        @test isempty(trim_angle(solver, body_aero, beta, wind_speed;
+                                 alpha_range=deg2rad.(4:2:12)))
+    end
+end