diff --git a/benchmarks/benchmark_against_fortran_run.jl b/benchmarks/benchmark_against_fortran_run.jl
index 8d0fb078..b48cb39c 100644
--- a/benchmarks/benchmark_against_fortran_run.jl
+++ b/benchmarks/benchmark_against_fortran_run.jl
@@ -55,9 +55,9 @@ import GeneralizedPerturbedEquilibrium.Analysis.PerturbedEquilibriumModes
 # ─── Argument parsing ────────────────────────────────────────────────────────
 
 function parse_args(args)
-    fortran_dir     = ""
-    do_plot         = false
-    skip_run        = false
+    fortran_dir = ""
+    do_plot = false
+    skip_run = false
     use_file_forcing = false
     i = 1
     while i <= length(args)
@@ -78,7 +78,7 @@ function parse_args(args)
     isempty(fortran_dir) && error(
         "Usage: benchmark_fortran.jl <fortran_dir> [--plot] [--skip-run] [--use-file-forcing]")
     return (fortran_dir=abspath(fortran_dir), do_plot=do_plot,
-            skip_run=skip_run, use_file_forcing=use_file_forcing)
+        skip_run=skip_run, use_file_forcing=use_file_forcing)
 end
 
 # ─── Fortran namelist parsing ─────────────────────────────────────────────────
@@ -109,7 +109,7 @@ end
 
 function _find_indexed2_float(text::String, name::String, i::Int, j::Int; default=nothing)
     m = match(Regex(
-        "\\b$(name)\\s*\\(\\s*$(i)\\s*,\\s*$(j)\\s*\\)\\s*=\\s*([+-]?[\\d.]+(?:[eE][+-]?\\d+)?)", "i"),
+            "\\b$(name)\\s*\\(\\s*$(i)\\s*,\\s*$(j)\\s*\\)\\s*=\\s*([+-]?[\\d.]+(?:[eE][+-]?\\d+)?)", "i"),
         text)
     isnothing(m) ? default : parse(Float64, m.captures[1])
 end
@@ -140,34 +140,34 @@ struct FortranRunParams
 end
 
 function parse_fortran_run(dir::String)
-    eq_text   = _strip_fortran_comments(read(joinpath(dir, "equil.in"), String))
-    dcon_text = _strip_fortran_comments(read(joinpath(dir, "dcon.in"),  String))
+    eq_text = _strip_fortran_comments(read(joinpath(dir, "equil.in"), String))
+    dcon_text = _strip_fortran_comments(read(joinpath(dir, "dcon.in"), String))
     coil_text = isfile(joinpath(dir, "coil.in")) ?
-        _strip_fortran_comments(read(joinpath(dir, "coil.in"), String)) : ""
+                _strip_fortran_comments(read(joinpath(dir, "coil.in"), String)) : ""
 
-    eq_type   = _find_string(eq_text, "eq_type";   default="efit")
-    eq_file   = _find_string(eq_text, "eq_filename"; default="")
-    jac_type  = _find_string(eq_text, "jac_type";  default="hamada")
+    eq_type = _find_string(eq_text, "eq_type"; default="efit")
+    eq_file = _find_string(eq_text, "eq_filename"; default="")
+    jac_type = _find_string(eq_text, "jac_type"; default="hamada")
     grid_type = _find_string(eq_text, "grid_type"; default="ldp")
-    psilow    = _find_scalar(eq_text, "psilow";    default=1e-4)
-    psihigh   = _find_scalar(eq_text, "psihigh";   default=0.993)
-    mpsi      = _find_int(eq_text,    "mpsi";      default=128)
-    mtheta_eq = _find_int(eq_text,    "mtheta";    default=256)
+    psilow = _find_scalar(eq_text, "psilow"; default=1e-4)
+    psihigh = _find_scalar(eq_text, "psihigh"; default=0.993)
+    mpsi = _find_int(eq_text, "mpsi"; default=128)
+    mtheta_eq = _find_int(eq_text, "mtheta"; default=256)
 
-    nn          = _find_int(dcon_text, "nn";           default=1)
-    delta_mlow  = _find_int(dcon_text, "delta_mlow";  default=8)
+    nn = _find_int(dcon_text, "nn"; default=1)
+    delta_mlow = _find_int(dcon_text, "delta_mlow"; default=8)
     delta_mhigh = _find_int(dcon_text, "delta_mhigh"; default=8)
-    dmlim       = _find_scalar(dcon_text, "dmlim";    default=0.2)
-    qlow        = _find_scalar(dcon_text, "qlow";     default=1.02)
-    psiedge     = _find_scalar(dcon_text, "psiedge";  default=1.0)
+    dmlim = _find_scalar(dcon_text, "dmlim"; default=0.2)
+    qlow = _find_scalar(dcon_text, "qlow"; default=1.02)
+    psiedge = _find_scalar(dcon_text, "psiedge"; default=1.0)
     euler_stride = _find_int(dcon_text, "euler_stride"; default=1)
 
-    mtheta_coil = _find_int(coil_text, "cmtheta";  default=480)
-    nzeta_coil  = _find_int(coil_text, "cmzeta";   default=40)
-    machine     = _find_string(coil_text, "machine"; default="d3d")
-    coil_num    = _find_int(coil_text, "coil_num"; default=0)
+    mtheta_coil = _find_int(coil_text, "cmtheta"; default=480)
+    nzeta_coil = _find_int(coil_text, "cmzeta"; default=40)
+    machine = _find_string(coil_text, "machine"; default="d3d")
+    coil_num = _find_int(coil_text, "coil_num"; default=0)
 
-    coil_names    = String[]
+    coil_names = String[]
     coil_currents = Vector{Float64}[]
     for ci in 1:coil_num
         name = _find_indexed_string(coil_text, "coil_name", ci; default="coil$ci")
@@ -189,11 +189,11 @@ end
 # ─── Load Fortran NC outputs ─────────────────────────────────────────────────
 
 function load_fortran_outputs(fortran_dir::String, nn::Int)
-    dcon_file    = joinpath(fortran_dir, "dcon_output_n$(nn).nc")
+    dcon_file = joinpath(fortran_dir, "dcon_output_n$(nn).nc")
     profile_file = joinpath(fortran_dir, "gpec_profile_output_n$(nn).nc")
     control_file = joinpath(fortran_dir, "gpec_control_output_n$(nn).nc")
 
-    isfile(dcon_file)    || error("Not found: $dcon_file")
+    isfile(dcon_file) || error("Not found: $dcon_file")
     isfile(profile_file) || error("Not found: $profile_file")
     isfile(control_file) || error("Not found: $control_file")
 
@@ -201,16 +201,16 @@ function load_fortran_outputs(fortran_dir::String, nn::Int)
 
     NCDataset(dcon_file, "r") do ds
         fort["psilim"] = Float64(ds.attrib["psilim"])
-        fort["qlim"]   = Float64(ds.attrib["qlim"])
-        fort["mlow"]   = Int(ds.attrib["mlow"])
-        fort["mhigh"]  = Int(ds.attrib["mhigh"])
-        fort["mpert"]  = Int(ds.attrib["mpert"])
+        fort["qlim"] = Float64(ds.attrib["qlim"])
+        fort["mlow"] = Int(ds.attrib["mlow"])
+        fort["mhigh"] = Int(ds.attrib["mhigh"])
+        fort["mpert"] = Int(ds.attrib["mpert"])
 
         fort["psi_n_dcon"] = Float64.(ds["psi_n"][:])
-        fort["q_dcon"]     = Float64.(ds["q"][:])
-        fort["di"]         = Float64.(ds["di"][:])
-        fort["dr"]         = Float64.(ds["dr"][:])
-        fort["m_vals"]     = Int.(ds["m"][:])
+        fort["q_dcon"] = Float64.(ds["q"][:])
+        fort["di"] = Float64.(ds["di"][:])
+        fort["dr"] = Float64.(ds["dr"][:])
+        fort["m_vals"] = Int.(ds["m"][:])
 
         wt_ev = Float64.(ds["W_t_eigenvalue"][:, :])
         fort["W_t_eigenvalue"] = complex.(wt_ev[:, 1], wt_ev[:, 2])
@@ -218,7 +218,7 @@ function load_fortran_outputs(fortran_dir::String, nn::Int)
 
     NCDataset(profile_file, "r") do ds
         fort["psi_n_rational"] = Float64.(ds["psi_n_rational"][:])
-        fort["q_rational"]     = Float64.(ds["q_rational"][:])
+        fort["q_rational"] = Float64.(ds["q_rational"][:])
 
         phi_res = Float64.(ds["Phi_res"][:, :])
         fort["Phi_res"] = complex.(phi_res[:, 1], phi_res[:, 2])
@@ -237,7 +237,7 @@ function load_fortran_outputs(fortran_dir::String, nn::Int)
 
         # radial profiles vs psi for all modes
         fort["psi_n_prof"] = Float64.(ds["psi_n"][:])
-        fort["m_out"]      = Int.(ds["m_out"][:])
+        fort["m_out"] = Int.(ds["m_out"][:])
         b_n_raw = Float64.(ds["b_n"][:, :, :])
         fort["b_n"] = complex.(b_n_raw[:, :, 1], b_n_raw[:, :, 2])
         jbgp_raw = Float64.(ds["Jbgradpsi"][:, :, :])
@@ -265,13 +265,13 @@ function load_fortran_outputs(fortran_dir::String, nn::Int)
     NCDataset(control_file, "r") do ds
         # Phi_coil: NC dim (i=2, coil_index=ncoil, m=mpert)
         phi_coil = Float64.(ds["Phi_coil"][:, :, :])
-        amp_real = dropdims(sum(phi_coil[:, :, 1], dims=2), dims=2)
-        amp_imag = dropdims(sum(phi_coil[:, :, 2], dims=2), dims=2)
+        amp_real = dropdims(sum(phi_coil[:, :, 1]; dims=2); dims=2)
+        amp_imag = dropdims(sum(phi_coil[:, :, 2]; dims=2); dims=2)
         fort["Phi_coil_amp"] = complex.(amp_real, amp_imag)
 
         phi_x = Float64.(ds["Phi_x"][:, :])
         fort["Phi_x_amp"] = abs.(complex.(phi_x[:, 1], phi_x[:, 2]))
-        fort["Phi_x"]     = complex.(phi_x[:, 1], phi_x[:, 2])
+        fort["Phi_x"] = complex.(phi_x[:, 1], phi_x[:, 2])
 
         phi_tot = Float64.(ds["Phi"][:, :])
         fort["Phi_tot"] = complex.(phi_tot[:, 1], phi_tot[:, 2])
@@ -295,7 +295,7 @@ function load_fortran_outputs(fortran_dir::String, nn::Int)
         if haskey(ds, "xi_r")
             ntheta = length(ds["theta_dcon"]) - 1  # drop duplicate endpoint
             for (key, nc_var) in [("xi_r_fun", "xi_r"), ("xi_z_fun", "xi_z"), ("xi_phi_fun", "xi_phi"),
-                                  ("b_r_fun",  "b_r"),  ("b_z_fun",  "b_z"),  ("b_phi_fun",  "b_phi")]
+                ("b_r_fun", "b_r"), ("b_z_fun", "b_z"), ("b_phi_fun", "b_phi")]
                 raw = Float64.(ds[nc_var][:, :, :])  # (psi, theta+1, 2)
                 fort[key] = complex.(raw[:, 1:ntheta, 1], raw[:, 1:ntheta, 2])
             end
@@ -313,50 +313,50 @@ function load_julia_outputs(h5_path::String)
     julia = Dict{String,Any}()
     h5open(h5_path, "r") do f
         julia["psilim"] = read(f, "info/psilim")
-        julia["qlim"]   = read(f, "info/qlim")
-        julia["et"]     = read(f, "vacuum/et")
-        julia["psi_q"]  = read(f, "splines/profiles/xs")
-        julia["q"]      = read(f, "splines/profiles/q")
-        julia["di"]     = haskey(f, "locstab/di") ? read(f, "locstab/di") : Float64[]
-        julia["dr"]     = haskey(f, "locstab/dr") ? read(f, "locstab/dr") : Float64[]
+        julia["qlim"] = read(f, "info/qlim")
+        julia["et"] = read(f, "vacuum/et")
+        julia["psi_q"] = read(f, "splines/profiles/xs")
+        julia["q"] = read(f, "splines/profiles/q")
+        julia["di"] = haskey(f, "locstab/di") ? read(f, "locstab/di") : Float64[]
+        julia["dr"] = haskey(f, "locstab/dr") ? read(f, "locstab/dr") : Float64[]
 
         julia["psio"] = haskey(f, "equil/psio") ? read(f, "equil/psio") : NaN
 
         sc = "perturbed_equilibrium/singular_coupling"
-        julia["rational_psi"]       = haskey(f, "$sc/rational_psi")       ? read(f, "$sc/rational_psi")       : Float64[]
-        julia["rational_q"]         = haskey(f, "$sc/rational_q")         ? read(f, "$sc/rational_q")         : Float64[]
-        julia["rational_n"]         = haskey(f, "$sc/rational_n")         ? read(f, "$sc/rational_n")         : Int[]
-        julia["rational_m_res"]     = haskey(f, "$sc/rational_m_res")     ? read(f, "$sc/rational_m_res")     : Int[]
-        julia["resonant_area_weighted_field"]     = haskey(f, "$sc/resonant_area_weighted_field")     ? read(f, "$sc/resonant_area_weighted_field")     : ComplexF64[]
-        julia["resonant_current"]   = haskey(f, "$sc/resonant_current")   ? read(f, "$sc/resonant_current")   : ComplexF64[]
-        julia["island_half_width"]  = haskey(f, "$sc/island_half_width")  ? read(f, "$sc/island_half_width")  : Float64[]
+        julia["rational_psi"] = haskey(f, "$sc/rational_psi") ? read(f, "$sc/rational_psi") : Float64[]
+        julia["rational_q"] = haskey(f, "$sc/rational_q") ? read(f, "$sc/rational_q") : Float64[]
+        julia["rational_n"] = haskey(f, "$sc/rational_n") ? read(f, "$sc/rational_n") : Int[]
+        julia["rational_m_res"] = haskey(f, "$sc/rational_m_res") ? read(f, "$sc/rational_m_res") : Int[]
+        julia["resonant_area_weighted_field"] = haskey(f, "$sc/resonant_area_weighted_field") ? read(f, "$sc/resonant_area_weighted_field") : ComplexF64[]
+        julia["resonant_current"] = haskey(f, "$sc/resonant_current") ? read(f, "$sc/resonant_current") : ComplexF64[]
+        julia["island_half_width"] = haskey(f, "$sc/island_half_width") ? read(f, "$sc/island_half_width") : Float64[]
         julia["chirikov_parameter"] = haskey(f, "$sc/chirikov_parameter") ? read(f, "$sc/chirikov_parameter") : Float64[]
-        julia["delta_prime"]        = haskey(f, "$sc/delta_prime")        ? read(f, "$sc/delta_prime")        : ComplexF64[]
+        julia["delta_prime"] = haskey(f, "$sc/delta_prime") ? read(f, "$sc/delta_prime") : ComplexF64[]
 
         pe = "perturbed_equilibrium"
         # Fortran Phi_x/Phi_tot are the area-weighted field b̄ (tesla), matching forcing/response_b_area directly.
-        julia["forcing_vec"]  = haskey(f, "$pe/forcing_b_area")  ? read(f, "$pe/forcing_b_area")  : ComplexF64[]   # Phi_x   = b̄_x
+        julia["forcing_vec"] = haskey(f, "$pe/forcing_b_area") ? read(f, "$pe/forcing_b_area") : ComplexF64[]   # Phi_x   = b̄_x
         julia["response_vec"] = haskey(f, "$pe/response_b_area") ? read(f, "$pe/response_b_area") : ComplexF64[]   # Phi_tot = b̄_tot
-        julia["b_n"]          = haskey(f, "$pe/response/b_n")        ? read(f, "$pe/response/b_n")        : Matrix{ComplexF64}(undef, 0, 0)
-        julia["Jbgradpsi"]    = haskey(f, "$pe/response/b_psi_area_weighted") ? read(f, "$pe/response/b_psi_area_weighted") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["xi_psi"]       = haskey(f, "$pe/response/xi_psi")     ? read(f, "$pe/response/xi_psi")     : Matrix{ComplexF64}(undef, 0, 0)
-        julia["xi_n"]         = haskey(f, "$pe/response/xi_n")       ? read(f, "$pe/response/xi_n")       : Matrix{ComplexF64}(undef, 0, 0)
-        julia["clebsch_psi1"] = haskey(f, "$pe/response/clebsch_psi1")  ? read(f, "$pe/response/clebsch_psi1")  : Matrix{ComplexF64}(undef, 0, 0)
-        julia["clebsch_alpha"]= haskey(f, "$pe/response/clebsch_alpha") ? read(f, "$pe/response/clebsch_alpha") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["psi_grid"]     = haskey(f, "integration/psi")    ? read(f, "integration/psi")    : Float64[]
+        julia["b_n"] = haskey(f, "$pe/response/b_n") ? read(f, "$pe/response/b_n") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["Jbgradpsi"] = haskey(f, "$pe/response/b_psi_area_weighted") ? read(f, "$pe/response/b_psi_area_weighted") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["xi_psi"] = haskey(f, "$pe/response/xi_psi") ? read(f, "$pe/response/xi_psi") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["xi_n"] = haskey(f, "$pe/response/xi_n") ? read(f, "$pe/response/xi_n") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["clebsch_psi1"] = haskey(f, "$pe/response/clebsch_psi1") ? read(f, "$pe/response/clebsch_psi1") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["clebsch_alpha"] = haskey(f, "$pe/response/clebsch_alpha") ? read(f, "$pe/response/clebsch_alpha") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["psi_grid"] = haskey(f, "integration/psi") ? read(f, "integration/psi") : Float64[]
 
         # R,Z,φ: loaded via modes_to_theta helper below (not raw modes)
         julia["h5_path"] = h5_path  # stash for modes_to_theta
         # mn_index[:, 1] = m values, mn_index[:, 2] = n values for each mode index
-        julia["m_modes"]      = haskey(f, "info/mn_index") ? Int.(read(f, "info/mn_index")[:, 1]) : Int[]
+        julia["m_modes"] = haskey(f, "info/mn_index") ? Int.(read(f, "info/mn_index")[:, 1]) : Int[]
 
         # Control surface matrices
         rm = "$pe/response_matrices"
-        julia["permeability"]       = haskey(f, "$rm/permeability")       ? read(f, "$rm/permeability")       : Matrix{ComplexF64}(undef, 0, 0)
-        julia["plasma_inductance"]  = haskey(f, "$rm/plasma_inductance")  ? read(f, "$rm/plasma_inductance")  : Matrix{ComplexF64}(undef, 0, 0)
+        julia["permeability"] = haskey(f, "$rm/permeability") ? read(f, "$rm/permeability") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["plasma_inductance"] = haskey(f, "$rm/plasma_inductance") ? read(f, "$rm/plasma_inductance") : Matrix{ComplexF64}(undef, 0, 0)
         julia["surface_inductance"] = haskey(f, "$rm/surface_inductance") ? read(f, "$rm/surface_inductance") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["reluctance"]         = haskey(f, "$rm/reluctance")         ? read(f, "$rm/reluctance")         : Matrix{ComplexF64}(undef, 0, 0)
-        julia["wt0"]                = haskey(f, "vacuum/wt0")             ? read(f, "vacuum/wt0")             : Matrix{ComplexF64}(undef, 0, 0)
+        julia["reluctance"] = haskey(f, "$rm/reluctance") ? read(f, "$rm/reluctance") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["wt0"] = haskey(f, "vacuum/wt0") ? read(f, "vacuum/wt0") : Matrix{ComplexF64}(undef, 0, 0)
     end
     return julia
 end
@@ -410,8 +410,6 @@ function write_gpec_toml_coil(
         println(io)
         println(io, "[ForceFreeStates]")
         println(io, "local_stability_flag = true")
-        println(io, "mat_flag = true")
-        println(io, "ode_flag = true")
         println(io, "vac_flag = true")
         @printf(io, "qlow     = %.4f\n", p.qlow)
         # Match the Fortran sas_flag truncation: integration stops at q = qhigh
@@ -475,8 +473,6 @@ function write_gpec_toml_file(
         println(io)
         println(io, "[ForceFreeStates]")
         println(io, "local_stability_flag = true")
-        println(io, "mat_flag = true")
-        println(io, "ode_flag = true")
         println(io, "vac_flag = true")
         @printf(io, "qlow     = %.4f\n", p.qlow)
         @printf(io, "psiedge  = %.4f\n", p.psiedge)
@@ -506,15 +502,15 @@ end
 
 function setup_equil(p::FortranRunParams, fortran_dir::String)
     eq_dict = Dict{String,Any}(
-        "eq_type"      => p.eq_type,
-        "eq_filename"  => p.eq_filename,
-        "jac_type"     => p.jac_type,
-        "psilow"       => p.psilow,
-        "psihigh"      => p.psihigh,
-        "mpsi"         => p.mpsi,
-        "mtheta"       => p.mtheta_equil,
-        "grid_type"    => p.grid_type,
-        "etol"         => 1e-7,
+        "eq_type" => p.eq_type,
+        "eq_filename" => p.eq_filename,
+        "jac_type" => p.jac_type,
+        "psilow" => p.psilow,
+        "psihigh" => p.psihigh,
+        "mpsi" => p.mpsi,
+        "mtheta" => p.mtheta_equil,
+        "grid_type" => p.grid_type,
+        "etol" => 1e-7
     )
     eq_config = Equilibrium.EquilibriumConfig(eq_dict, fortran_dir)
     return Equilibrium.setup_equilibrium(eq_config)
@@ -537,8 +533,8 @@ function compare_phix(
 
     psilim = fort["psilim"]
     m_vals = fort["m_vals"]
-    mlow   = minimum(m_vals)
-    mhigh  = maximum(m_vals)
+    mlow = minimum(m_vals)
+    mhigh = maximum(m_vals)
 
     coil_sets = ForcingTerms.CoilSet[]
     for (name, currents) in zip(p.coil_names, p.coil_currents)
@@ -551,21 +547,25 @@ function compare_phix(
 
     grid = ForcingTerms.sample_boundary_grid(equil, p.mtheta_coil, p.nzeta_coil; psi=psilim)
     nobs = p.mtheta_coil * p.nzeta_coil
-    obs_R   = zeros(nobs); obs_phi = zeros(nobs); obs_Z = zeros(nobs)
+    obs_R = zeros(nobs);
+    obs_phi = zeros(nobs);
+    obs_Z = zeros(nobs)
     for j in 1:p.nzeta_coil, i in 1:p.mtheta_coil
         idx = i + (j - 1) * p.mtheta_coil
-        obs_R[idx]   = grid.R[i]
+        obs_R[idx] = grid.R[i]
         obs_phi[idx] = grid.phi_grid[j] + grid.phi_offset[i]
-        obs_Z[idx]   = grid.Z[i]
+        obs_Z[idx] = grid.Z[i]
     end
-    B_R = zeros(nobs); B_phi = zeros(nobs); B_Z = zeros(nobs)
+    B_R = zeros(nobs);
+    B_phi = zeros(nobs);
+    B_Z = zeros(nobs)
     ForcingTerms.compute_biot_savart_boundary!(B_R, B_phi, B_Z, obs_R, obs_phi, obs_Z, coil_sets)
     bn = zeros(p.mtheta_coil, p.nzeta_coil)
     ForcingTerms.project_normal_flux!(bn, B_R, B_Z, grid)
     julia_modes = ForcingTerms.fourier_decompose_bn(bn, grid, p.nn, mlow, mhigh)
 
     julia_amps = Dict(md.m => abs(md.amplitude) for md in julia_modes)
-    fort_amps  = Dict(m => fort["Phi_x_amp"][k] for (k, m) in enumerate(m_vals))
+    fort_amps = Dict(m => fort["Phi_x_amp"][k] for (k, m) in enumerate(m_vals))
 
     return julia_amps, fort_amps, m_vals, equil
 end
@@ -575,7 +575,7 @@ end
 rel_pct(a, b) = abs(b) > 1e-30 ? @sprintf("%.1f%%", 100 * abs(a - b) / abs(b)) : "N/A"
 
 function step_series(m_vals, amps)
-    m_ext   = [m_vals[1] - 1; m_vals; m_vals[end] + 1]
+    m_ext = [m_vals[1] - 1; m_vals; m_vals[end] + 1]
     amp_ext = [0.0; amps; 0.0]
     return m_ext, amp_ext
 end
@@ -605,7 +605,7 @@ function build_comparison_table(fort, julia, fortran_dir, bench_dir, nn)
     push!(lines, @sprintf("%-12s  %10s  %10s  %10s  %6s", "", "Fortran", "Julia", "AbsDiff", "Rel%"))
     for (label, fv, jv) in [
         ("psilim", fort["psilim"], julia["psilim"]),
-        ("qlim",   fort["qlim"],   julia["qlim"]),
+        ("qlim", fort["qlim"], julia["qlim"])
     ]
         push!(lines, @sprintf("%-12s  %10.6f  %10.6f  %10.6f  %6s",
             label, fv, jv, abs(fv - jv), rel_pct(fv, jv)))
@@ -624,12 +624,12 @@ function build_comparison_table(fort, julia, fortran_dir, bench_dir, nn)
 
     push!(lines, "--- q Profile (RMS on common grid) ---")
     fq_psi = fort["psi_n_dcon"]
-    fq     = fort["q_dcon"]
+    fq = fort["q_dcon"]
     jq_psi = julia["psi_q"]
-    jq     = julia["q"]
+    jq = julia["q"]
     if !isempty(jq) && !isempty(fq)
         jq_on_fgrid = [_interp1(jq_psi, jq, p) for p in fq_psi]
-        rms = sqrt(mean((fq .- jq_on_fgrid).^2))
+        rms = sqrt(mean((fq .- jq_on_fgrid) .^ 2))
         push!(lines, @sprintf("RMS(q_fort - q_julia) = %.5f  (%.2f%% relative)",
             rms, rms / mean(abs.(fq)) * 100))
     else
@@ -638,42 +638,47 @@ function build_comparison_table(fort, julia, fortran_dir, bench_dir, nn)
     push!(lines, "")
 
     push!(lines, "--- Mercier Criterion ---")
-    fdi = fort["di"]; fdr = fort["dr"]
-    jdi = julia["di"]; jdr = julia["dr"]
+    fdi = fort["di"];
+    fdr = fort["dr"]
+    jdi = julia["di"];
+    jdr = julia["dr"]
     if !isempty(jdi) && !isempty(fdi)
         jdi_on_fg = [_interp1(jq_psi, jdi, p) for p in fq_psi]
         jdr_on_fg = !isempty(jdr) ?
-            [_interp1(jq_psi, jdr, p) for p in fq_psi] : fill(NaN, length(fq_psi))
-        push!(lines, @sprintf("RMS(di_fort - di_julia) = %.5e", sqrt(mean((fdi .- jdi_on_fg).^2))))
-        push!(lines, @sprintf("RMS(dr_fort - dr_julia) = %.5e", sqrt(mean((fdr .- jdr_on_fg).^2))))
+                    [_interp1(jq_psi, jdr, p) for p in fq_psi] : fill(NaN, length(fq_psi))
+        push!(lines, @sprintf("RMS(di_fort - di_julia) = %.5e", sqrt(mean((fdi .- jdi_on_fg) .^ 2))))
+        push!(lines, @sprintf("RMS(dr_fort - dr_julia) = %.5e", sqrt(mean((fdr .- jdr_on_fg) .^ 2))))
     else
         push!(lines, "Mercier data unavailable")
     end
     push!(lines, "")
 
     push!(lines, "--- Phi_x and Phi_tot Spectrum (n=$nn) ---")
-    m_ctrl   = fort["m_vals_ctrl"]
-    f_phi_x  = fort["Phi_x"]
+    m_ctrl = fort["m_vals_ctrl"]
+    f_phi_x = fort["Phi_x"]
     f_phi_tot = fort["Phi_tot"]
-    j_fvec   = julia["forcing_vec"]
-    j_rvec   = julia["response_vec"]
+    j_fvec = julia["forcing_vec"]
+    j_rvec = julia["response_vec"]
     j_mmodes = julia["m_modes"]
     if !isempty(j_fvec) && !isempty(f_phi_x)
         push!(lines, @sprintf("%-5s  %-14s  %-14s  %-6s  %-14s  %-14s  %-6s",
             "m", "|Phi_x_fort|", "|Phi_x_julia|", "Rel%", "|Phi_tot_fort|", "|Phi_tot_julia|", "Rel%"))
         for (k, m) in enumerate(m_ctrl)
             j_idx = findfirst(==(m), j_mmodes)
-            j_fx  = isnothing(j_idx) ? NaN : abs(j_fvec[j_idx])
-            j_ft  = isnothing(j_idx) ? NaN : abs(j_rvec[j_idx])
-            f_fx  = abs(f_phi_x[k])
-            f_ft  = abs(f_phi_tot[k])
-            push!(lines, @sprintf("%-5d  %14.4e  %14s  %-6s  %14.4e  %14s  %-6s",
-                m, f_fx,
-                isnan(j_fx) ? "N/A" : @sprintf("%.4e", j_fx),
-                isnan(j_fx) ? "N/A" : rel_pct(f_fx, j_fx),
-                f_ft,
-                isnan(j_ft) ? "N/A" : @sprintf("%.4e", j_ft),
-                isnan(j_ft) ? "N/A" : rel_pct(f_ft, j_ft)))
+            j_fx = isnothing(j_idx) ? NaN : abs(j_fvec[j_idx])
+            j_ft = isnothing(j_idx) ? NaN : abs(j_rvec[j_idx])
+            f_fx = abs(f_phi_x[k])
+            f_ft = abs(f_phi_tot[k])
+            push!(
+                lines,
+                @sprintf("%-5d  %14.4e  %14s  %-6s  %14.4e  %14s  %-6s",
+                    m, f_fx,
+                    isnan(j_fx) ? "N/A" : @sprintf("%.4e", j_fx),
+                    isnan(j_fx) ? "N/A" : rel_pct(f_fx, j_fx),
+                    f_ft,
+                    isnan(j_ft) ? "N/A" : @sprintf("%.4e", j_ft),
+                    isnan(j_ft) ? "N/A" : rel_pct(f_ft, j_ft))
+            )
         end
     else
         push!(lines, "Phi_tot data unavailable")
@@ -737,40 +742,46 @@ function build_comparison_table(fort, julia, fortran_dir, bench_dir, nn)
     push!(lines, "")
 
     push!(lines, "--- Resonant Surface Comparison (n=$nn) ---")
-    push!(lines, @sprintf(
-        "%-4s  %-8s  %-8s  %-14s  %-14s  %-6s  %-10s  %-10s  %-6s  %-7s  %-7s  %-6s",
-        "q", "psi_fort", "psi_julia", "|Phi_res_fort|", "|Phi_res_julia|", "Rel%",
-        "w_isl_fort", "w_isl_julia", "Rel%", "K_fort", "K_julia", "Rel%"))
+    push!(
+        lines,
+        @sprintf(
+            "%-4s  %-8s  %-8s  %-14s  %-14s  %-6s  %-10s  %-10s  %-6s  %-7s  %-7s  %-6s",
+            "q", "psi_fort", "psi_julia", "|Phi_res_fort|", "|Phi_res_julia|", "Rel%",
+            "w_isl_fort", "w_isl_julia", "Rel%", "K_fort", "K_julia", "Rel%")
+    )
 
     f_psi_rat = fort["psi_n_rational"]
-    f_q_rat   = fort["q_rational"]
+    f_q_rat = fort["q_rational"]
     f_phi_res = fort["Phi_res"]
-    f_w_isl   = fort["w_isl"]
-    f_K_isl   = fort["K_isl"]
+    f_w_isl = fort["w_isl"]
+    f_K_isl = fort["K_isl"]
 
     for (i, (fq_r, fp_r)) in enumerate(zip(f_q_rat, f_psi_rat))
         q_int = round(Int, fq_r)
-        row   = find_julia_row(julia, q_int, nn)
-
-        jp_r  = (row > 0 && !isempty(julia["rational_psi"]))     ? julia["rational_psi"][row]        : NaN
-        j_phi = (row > 0 && !isempty(julia["resonant_area_weighted_field"]))   ? abs(julia["resonant_area_weighted_field"][row])  : NaN
-        j_wisl= (row > 0 && !isempty(julia["island_half_width"])) ? 2*julia["island_half_width"][row] : NaN
-        j_kch = (row > 0 && !isempty(julia["chirikov_parameter"])) ? julia["chirikov_parameter"][row]  : NaN
-
-        push!(lines, @sprintf(
-            "%-4d  %8.4f  %8s  %14.4e  %14s  %-6s  %10.4e  %10s  %-6s  %7.4f  %7s  %-6s",
-            q_int, fp_r,
-            isnan(jp_r)   ? "N/A" : @sprintf("%.4f", jp_r),
-            abs(f_phi_res[i]),
-            isnan(j_phi)  ? "N/A" : @sprintf("%.4e", j_phi),
-            isnan(j_phi)  ? "N/A" : rel_pct(abs(f_phi_res[i]), j_phi),
-            f_w_isl[i],
-            isnan(j_wisl) ? "N/A" : @sprintf("%.4e", j_wisl),
-            isnan(j_wisl) ? "N/A" : rel_pct(f_w_isl[i], j_wisl),
-            f_K_isl[i],
-            isnan(j_kch)  ? "N/A" : @sprintf("%.4f", j_kch),
-            isnan(j_kch)  ? "N/A" : rel_pct(f_K_isl[i], j_kch),
-        ))
+        row = find_julia_row(julia, q_int, nn)
+
+        jp_r = (row > 0 && !isempty(julia["rational_psi"])) ? julia["rational_psi"][row] : NaN
+        j_phi = (row > 0 && !isempty(julia["resonant_area_weighted_field"])) ? abs(julia["resonant_area_weighted_field"][row]) : NaN
+        j_wisl = (row > 0 && !isempty(julia["island_half_width"])) ? 2*julia["island_half_width"][row] : NaN
+        j_kch = (row > 0 && !isempty(julia["chirikov_parameter"])) ? julia["chirikov_parameter"][row] : NaN
+
+        push!(
+            lines,
+            @sprintf(
+                "%-4d  %8.4f  %8s  %14.4e  %14s  %-6s  %10.4e  %10s  %-6s  %7.4f  %7s  %-6s",
+                q_int, fp_r,
+                isnan(jp_r) ? "N/A" : @sprintf("%.4f", jp_r),
+                abs(f_phi_res[i]),
+                isnan(j_phi) ? "N/A" : @sprintf("%.4e", j_phi),
+                isnan(j_phi) ? "N/A" : rel_pct(abs(f_phi_res[i]), j_phi),
+                f_w_isl[i],
+                isnan(j_wisl) ? "N/A" : @sprintf("%.4e", j_wisl),
+                isnan(j_wisl) ? "N/A" : rel_pct(f_w_isl[i], j_wisl),
+                f_K_isl[i],
+                isnan(j_kch) ? "N/A" : @sprintf("%.4f", j_kch),
+                isnan(j_kch) ? "N/A" : rel_pct(f_K_isl[i], j_kch),
+            )
+        )
     end
 
     push!(lines, "")
@@ -783,32 +794,36 @@ function build_comparison_table(fort, julia, fortran_dir, bench_dir, nn)
     push!(lines, @sprintf("%-4s  %-8s  %-14s  %-14s  %-8s",
         "q", "psi_fort", "j_c_fort", "j_c_julia", "ratio j/f"))
 
-    f_delta  = fort["Delta"]
-    f_i_res  = fort["I_res"]
-    j_delta  = julia["delta_prime"]
-    j_i_res  = julia["resonant_current"]
-    j_m_res  = julia["rational_m_res"]
+    f_delta = fort["Delta"]
+    f_i_res = fort["I_res"]
+    j_delta = julia["delta_prime"]
+    j_i_res = julia["resonant_current"]
+    j_m_res = julia["rational_m_res"]
 
     for (i, (fq_r, fp_r)) in enumerate(zip(f_q_rat, f_psi_rat))
-        q_int  = round(Int, fq_r)
+        q_int = round(Int, fq_r)
         m_res_fort = q_int * nn  # m_res = q * n at rational surface
         row = find_julia_row(julia, q_int, nn)
 
         jc_fort = (i <= length(f_delta) && i <= length(f_i_res) &&
                    abs(f_delta[i]) > 1e-30 && chi1 > 0) ?
-            abs(f_i_res[i]) * abs(m_res_fort) / (abs(f_delta[i]) * chi1) : NaN
+                  abs(f_i_res[i]) * abs(m_res_fort) / (abs(f_delta[i]) * chi1) : NaN
 
-        jc_julia = (row > 0 && row <= length(j_delta) && row <= length(j_i_res) &&
-                    row <= length(j_m_res) && abs(j_delta[row]) > 1e-30 && chi1 > 0) ?
+        jc_julia =
+            (row > 0 && row <= length(j_delta) && row <= length(j_i_res) &&
+             row <= length(j_m_res) && abs(j_delta[row]) > 1e-30 && chi1 > 0) ?
             abs(j_i_res[row]) * abs(j_m_res[row]) / (abs(j_delta[row]) * chi1) : NaN
 
         ratio = (!isnan(jc_fort) && !isnan(jc_julia) && jc_fort > 0) ? jc_julia / jc_fort : NaN
 
-        push!(lines, @sprintf("%-4d  %8.4f  %14s  %14s  %8s",
-            q_int, fp_r,
-            isnan(jc_fort)  ? "N/A" : @sprintf("%.4e", jc_fort),
-            isnan(jc_julia) ? "N/A" : @sprintf("%.4e", jc_julia),
-            isnan(ratio)    ? "N/A" : @sprintf("%.4f", ratio)))
+        push!(
+            lines,
+            @sprintf("%-4d  %8.4f  %14s  %14s  %8s",
+                q_int, fp_r,
+                isnan(jc_fort) ? "N/A" : @sprintf("%.4e", jc_fort),
+                isnan(jc_julia) ? "N/A" : @sprintf("%.4e", jc_julia),
+                isnan(ratio) ? "N/A" : @sprintf("%.4f", ratio))
+        )
     end
 
     push!(lines, "")
@@ -818,19 +833,22 @@ function build_comparison_table(fort, julia, fortran_dir, bench_dir, nn)
         "|I_res_fort|", "|I_res_julia|", "Rel%"))
     for (i, (fq_r, fp_r)) in enumerate(zip(f_q_rat, f_psi_rat))
         q_int = round(Int, fq_r)
-        row   = find_julia_row(julia, q_int, nn)
+        row = find_julia_row(julia, q_int, nn)
         fd = (i <= length(f_delta)) ? abs(f_delta[i]) : NaN
         jd = (row > 0 && row <= length(j_delta)) ? abs(j_delta[row]) : NaN
         fi = (!isempty(f_i_res) && i <= length(f_i_res)) ? abs(f_i_res[i]) : NaN
         ji = (row > 0 && row <= length(j_i_res)) ? abs(j_i_res[row]) : NaN
-        push!(lines, @sprintf("%-4d  %8.4f  %14s  %14s  %8s  %14s  %14s  %8s",
-            q_int, fp_r,
-            isnan(fd) ? "N/A" : @sprintf("%.4e", fd),
-            isnan(jd) ? "N/A" : @sprintf("%.4e", jd),
-            (isnan(fd) || isnan(jd)) ? "N/A" : rel_pct(fd, jd),
-            isnan(fi) ? "N/A" : @sprintf("%.4e", fi),
-            isnan(ji) ? "N/A" : @sprintf("%.4e", ji),
-            (isnan(fi) || isnan(ji)) ? "N/A" : rel_pct(fi, ji)))
+        push!(
+            lines,
+            @sprintf("%-4d  %8.4f  %14s  %14s  %8s  %14s  %14s  %8s",
+                q_int, fp_r,
+                isnan(fd) ? "N/A" : @sprintf("%.4e", fd),
+                isnan(jd) ? "N/A" : @sprintf("%.4e", jd),
+                (isnan(fd) || isnan(jd)) ? "N/A" : rel_pct(fd, jd),
+                isnan(fi) ? "N/A" : @sprintf("%.4e", fi),
+                isnan(ji) ? "N/A" : @sprintf("%.4e", ji),
+                (isnan(fi) || isnan(ji)) ? "N/A" : rel_pct(fi, ji))
+        )
     end
 
     push!(lines, "=" ^ 72)
@@ -838,8 +856,9 @@ function build_comparison_table(fort, julia, fortran_dir, bench_dir, nn)
 end
 
 function _interp1(x::AbstractVector, y::AbstractVector, xi::Real)
-    n = length(x); n < 2 && return y[1]
-    xi <= x[1]   && return y[1]
+    n = length(x);
+    n < 2 && return y[1]
+    xi <= x[1] && return y[1]
     xi >= x[end] && return y[end]
     i = clamp(searchsortedlast(x, xi), 1, n-1)
     t = (xi - x[i]) / (x[i+1] - x[i])
@@ -849,38 +868,38 @@ end
 # ─── Plot generation ─────────────────────────────────────────────────────────
 
 function generate_plots(fort, julia, bench_dir, nn)
-    f_q_rat  = fort["q_rational"]
-    n_surf   = length(f_q_rat)
+    f_q_rat = fort["q_rational"]
+    n_surf = length(f_q_rat)
 
     # --- Row 1: q profile, Mercier, W_t scatter, |P| diagonal ---
-    p1 = plot(fort["psi_n_dcon"], fort["q_dcon"],
+    p1 = plot(fort["psi_n_dcon"], fort["q_dcon"];
         label="Fortran", xlabel="ψ_n", ylabel="q", title="q profile", lw=2)
-    plot!(p1, julia["psi_q"], julia["q"], label="Julia", linestyle=:dash, lw=2)
+    plot!(p1, julia["psi_q"], julia["q"]; label="Julia", linestyle=:dash, lw=2)
 
-    p2 = plot(fort["psi_n_dcon"], fort["di"],
+    p2 = plot(fort["psi_n_dcon"], fort["di"];
         label="Fortran", xlabel="ψ_n", ylabel="D_I", title="Mercier D_I", lw=2)
-    !isempty(julia["di"]) && plot!(p2, julia["psi_q"], julia["di"], label="Julia", linestyle=:dash, lw=2)
-    hline!(p2, [0.0], color=:black, linestyle=:dot, label="")
+    !isempty(julia["di"]) && plot!(p2, julia["psi_q"], julia["di"]; label="Julia", linestyle=:dash, lw=2)
+    hline!(p2, [0.0]; color=:black, linestyle=:dot, label="")
 
-    p3 = plot(fort["psi_n_dcon"], fort["dr"],
+    p3 = plot(fort["psi_n_dcon"], fort["dr"];
         label="Fortran", xlabel="ψ_n", ylabel="D_R", title="Mercier D_R", lw=2)
-    !isempty(julia["dr"]) && plot!(p3, julia["psi_q"], julia["dr"], label="Julia", linestyle=:dash, lw=2)
-    hline!(p3, [0.0], color=:black, linestyle=:dot, label="")
+    !isempty(julia["dr"]) && plot!(p3, julia["psi_q"], julia["dr"]; label="Julia", linestyle=:dash, lw=2)
+    hline!(p3, [0.0]; color=:black, linestyle=:dot, label="")
 
     # W_t eigenvalue scatter (log y), all eigenvalues vs index
     f_wt = fort["W_t_eigenvalue"]
     j_et = julia["et"]
     p4 = plot(; xlabel="mode index", ylabel="|W_t|", title="|W_t| eigenvalues",
-               yscale=:log10, legend=:topright)
+        yscale=:log10, legend=:topright)
     if !isempty(f_wt)
         f_sorted = sort(abs.(f_wt))
         scatter!(p4, 1:length(f_sorted), f_sorted; color=:steelblue, markersize=5,
-                 markershape=:circle, label="Fortran")
+            markershape=:circle, label="Fortran")
     end
     if !isempty(j_et)
         j_sorted = sort(abs.(j_et))
         scatter!(p4, 1:length(j_sorted), j_sorted; color=:orange, markersize=5,
-                 markershape=:circle, label="Julia", markerstrokewidth=0)
+            markershape=:circle, label="Julia", markerstrokewidth=0)
     end
 
     # |P| diagonal comparison
@@ -888,12 +907,12 @@ function generate_plots(fort, julia, bench_dir, nn)
     j_P = julia["permeability"]
     m_ctrl = fort["m_vals_ctrl"]
     p5 = plot(; xlabel="Poloidal mode m", ylabel="|P[i,i]|",
-               title="|P| diagonal  (n=$nn)", legend=:topright)
+        title="|P| diagonal  (n=$nn)", legend=:topright)
     if !isempty(f_P) && !isempty(j_P)
         mpert_p = min(size(f_P, 1), size(j_P, 1), length(m_ctrl))
         ms = m_ctrl[1:mpert_p]
-        fp_diag = [abs(f_P[i,i]) for i in 1:mpert_p]
-        jp_diag = [abs(j_P[i,i]) for i in 1:mpert_p]
+        fp_diag = [abs(f_P[i, i]) for i in 1:mpert_p]
+        jp_diag = [abs(j_P[i, i]) for i in 1:mpert_p]
         me_f, fe = step_series(ms, fp_diag)
         me_j, je = step_series(ms, jp_diag)
         plot!(p5, me_f, fe; seriestype=:steppre, lw=2, color=:steelblue, label="Fortran")
@@ -905,35 +924,35 @@ function generate_plots(fort, julia, bench_dir, nn)
     # --- Row 2: Phi_x/tot spectrum and singular coupling vs psi ---
 
     # Phi_x & Phi_tot spectrum
-    m_ctrl   = fort["m_vals_ctrl"]
-    f_phi_x  = fort["Phi_x"]
+    m_ctrl = fort["m_vals_ctrl"]
+    f_phi_x = fort["Phi_x"]
     f_phi_tot = fort["Phi_tot"]
-    j_fvec   = julia["forcing_vec"]
-    j_rvec   = julia["response_vec"]
+    j_fvec = julia["forcing_vec"]
+    j_rvec = julia["response_vec"]
     j_mmodes = julia["m_modes"]
     p6 = plot(; xlabel="m", ylabel="|b̄| [T]", title="Phi_x & Phi_tot spectrum (n=$nn)",
-              legend=:topright)
+        legend=:topright)
     if !isempty(j_fvec) && !isempty(f_phi_x)
         jfx_amps = [abs(get(Dict(zip(j_mmodes, j_fvec)), m, NaN+0im)) for m in m_ctrl]
         jft_amps = [abs(get(Dict(zip(j_mmodes, j_rvec)), m, NaN+0im)) for m in m_ctrl]
         ffx_amps = abs.(f_phi_x)
         fft_amps = abs.(f_phi_tot)
         me, jfe = step_series(m_ctrl, replace(jfx_amps, NaN=>0.0))
-        _,  jte = step_series(m_ctrl, replace(jft_amps, NaN=>0.0))
-        _,  ffe = step_series(m_ctrl, ffx_amps)
-        _,  fte = step_series(m_ctrl, fft_amps)
+        _, jte = step_series(m_ctrl, replace(jft_amps, NaN=>0.0))
+        _, ffe = step_series(m_ctrl, ffx_amps)
+        _, fte = step_series(m_ctrl, fft_amps)
         plot!(p6, me, ffe; seriestype=:steppre, lw=2, color=:steelblue, label="Φ_x Fortran")
         plot!(p6, me, jfe; seriestype=:steppre, lw=2, color=:steelblue, linestyle=:dash, label="Φ_x Julia")
-        plot!(p6, me, fte; seriestype=:steppre, lw=2, color=:orange,    label="Φ_tot Fortran")
-        plot!(p6, me, jte; seriestype=:steppre, lw=2, color=:orange,    linestyle=:dash, label="Φ_tot Julia")
+        plot!(p6, me, fte; seriestype=:steppre, lw=2, color=:orange, label="Φ_tot Fortran")
+        plot!(p6, me, jte; seriestype=:steppre, lw=2, color=:orange, linestyle=:dash, label="Φ_tot Julia")
     else
         annotate!(p6, 0.5, 0.5, text("data unavailable", 10, :center))
     end
 
     # Helper: singular coupling vs psi scatter with q-value labels
     function _singcoup_panel(ylabel_str, title_str,
-                              f_psi, f_vals, f_q_labels,
-                              j_psi, j_vals)
+        f_psi, f_vals, f_q_labels,
+        j_psi, j_vals)
         p = plot(; xlabel="ψ_n", ylabel=ylabel_str, title=title_str, legend=:topright)
         if !isempty(f_psi) && !isempty(f_vals)
             scatter!(p, f_psi, f_vals; color=:steelblue, markersize=7, markershape=:circle, label="Fortran")
@@ -944,7 +963,7 @@ function generate_plots(fort, julia, bench_dir, nn)
         end
         if !isempty(j_psi) && !isempty(j_vals)
             scatter!(p, j_psi, j_vals; color=:orange, markersize=5, markershape=:circle,
-                     label="Julia", markerstrokewidth=0)
+                label="Julia", markerstrokewidth=0)
             plot!(p, j_psi, j_vals; color=:orange, lw=1, linestyle=:dot, label="")
         end
         return p
@@ -956,35 +975,35 @@ function generate_plots(fort, julia, bench_dir, nn)
 
     # Delta (Δ')
     f_delta_vals = isempty(fort["Delta"]) ? Float64[] : real.(fort["Delta"])
-    j_delta_psi  = julia["rational_psi"]
+    j_delta_psi = julia["rational_psi"]
     j_delta_vals = isempty(julia["delta_prime"]) ? Float64[] : real.(julia["delta_prime"])
     p7 = _singcoup_panel("Re(Δ')", "Δ' vs ψ  (n=$nn)",
-                          f_psi_rat, f_delta_vals, f_q_labels,
-                          j_delta_psi, j_delta_vals)
+        f_psi_rat, f_delta_vals, f_q_labels,
+        j_delta_psi, j_delta_vals)
 
     # I_res (resonant current)
     f_ires_vals = isempty(fort["I_res"]) ? Float64[] : abs.(fort["I_res"])
-    j_ires_psi  = julia["rational_psi"]
+    j_ires_psi = julia["rational_psi"]
     j_ires_vals = isempty(julia["resonant_current"]) ? Float64[] : abs.(julia["resonant_current"])
     p8 = _singcoup_panel("|I_res| [A]", "|I_res| vs ψ  (n=$nn)",
-                          f_psi_rat, f_ires_vals, f_q_labels,
-                          j_ires_psi, j_ires_vals)
+        f_psi_rat, f_ires_vals, f_q_labels,
+        j_ires_psi, j_ires_vals)
 
     # Phi_res (resonant flux / area)
     f_phires_vals = isempty(fort["Phi_res"]) ? Float64[] : abs.(fort["Phi_res"])
-    j_phires_psi  = julia["rational_psi"]
+    j_phires_psi = julia["rational_psi"]
     j_phires_vals = isempty(julia["resonant_area_weighted_field"]) ? Float64[] : abs.(julia["resonant_area_weighted_field"])
     p9 = _singcoup_panel("|Φ_res| [T]", "|Φ_res| vs ψ  (n=$nn)",
-                          f_psi_rat, f_phires_vals, f_q_labels,
-                          j_phires_psi, j_phires_vals)
+        f_psi_rat, f_phires_vals, f_q_labels,
+        j_phires_psi, j_phires_vals)
 
     # w_isl (island full width)
     f_wisl_vals = fort["w_isl"]
-    j_wisl_psi  = julia["rational_psi"]
+    j_wisl_psi = julia["rational_psi"]
     j_wisl_vals = isempty(julia["island_half_width"]) ? Float64[] : 2 .* julia["island_half_width"]
     p10 = _singcoup_panel("w_isl [ψ_n]", "Island width vs ψ  (n=$nn)",
-                           f_psi_rat, f_wisl_vals, f_q_labels,
-                           j_wisl_psi, j_wisl_vals)
+        f_psi_rat, f_wisl_vals, f_q_labels,
+        j_wisl_psi, j_wisl_vals)
 
     # --- Row 3: b_psi profiles for resonant m harmonics ---
     # Identify resonant m values from rational surfaces
@@ -994,23 +1013,23 @@ function generate_plots(fort, julia, bench_dir, nn)
 
     # Helper: build a normalised radial profile panel for one mode
     function _profile_panel(ylabel_str, title_str,
-                             f_psi, f_prof_complex, f_label,
-                             j_psi, j_prof_complex, j_label)
+        f_psi, f_prof_complex, f_label,
+        j_psi, j_prof_complex, j_label)
         pb = plot(; xlabel="ψ_n", ylabel=ylabel_str, title=title_str, legend=:topright)
         f_prof = abs.(f_prof_complex)
         j_prof = abs.(j_prof_complex)
-        f_max  = isempty(f_prof) ? 0.0 : maximum(f_prof)
-        j_max  = isempty(j_prof) ? 0.0 : maximum(j_prof)
+        f_max = isempty(f_prof) ? 0.0 : maximum(f_prof)
+        j_max = isempty(j_prof) ? 0.0 : maximum(j_prof)
         f_max > 0 && plot!(pb, f_psi, f_prof ./ f_max;
             lw=1.5, color=:steelblue, label=f_label,
             markershape=:x, markersize=2, markerstrokewidth=1, markeralpha=0.6)
         rms_diff = NaN
         if j_max > 0 && f_max > 0 && !isempty(j_psi)
-            j_on_fg  = [_interp1(j_psi, j_prof, p) for p in f_psi]
+            j_on_fg = [_interp1(j_psi, j_prof, p) for p in f_psi]
             j_max_on = maximum(j_on_fg)
-            f_norm   = f_prof ./ f_max
+            f_norm = f_prof ./ f_max
             if j_max_on > 0
-                rms_diff = sqrt(mean((j_on_fg ./ j_max_on .- f_norm).^2))
+                rms_diff = sqrt(mean((j_on_fg ./ j_max_on .- f_norm) .^ 2))
             end
             plot!(pb, j_psi, j_prof ./ j_max;
                 lw=1.5, color=:orange, linestyle=:dash, label=j_label,
@@ -1022,7 +1041,7 @@ function generate_plots(fort, julia, bench_dir, nn)
 
     j_psi_grid = julia["psi_grid"]
     f_psi_prof = fort["psi_n_prof"]
-    f_m_out    = fort["m_out"]
+    f_m_out = fort["m_out"]
 
     # Row 3: Jbgradpsi profiles — direct apples-to-apples (both = chi1*singfac*2πi*ξ_ψ in mode space)
     jbgp_panels = []
@@ -1033,10 +1052,13 @@ function generate_plots(fort, julia, bench_dir, nn)
         j_midx = findfirst(==(m_res), j_mmodes)
         fp = (!isnothing(f_midx) && !isempty(f_jbgp)) ? f_jbgp[:, f_midx] : ComplexF64[]
         jp = (!isnothing(j_midx) && !isempty(j_jbgp)) ? j_jbgp[:, j_midx] : ComplexF64[]
-        push!(jbgp_panels, _profile_panel(
-            "|Jbgradpsi| (norm.)", "Jbgradpsi  m=$m_res (n=$nn)",
-            isempty(fp) ? Float64[] : f_psi_prof, fp, "Fortran",
-            isempty(jp) ? Float64[] : j_psi_grid, jp, "Julia"))
+        push!(
+            jbgp_panels,
+            _profile_panel(
+                "|Jbgradpsi| (norm.)", "Jbgradpsi  m=$m_res (n=$nn)",
+                isempty(fp) ? Float64[] : f_psi_prof, fp, "Fortran",
+                isempty(jp) ? Float64[] : j_psi_grid, jp, "Julia")
+        )
     end
     while length(jbgp_panels) < 4
         push!(jbgp_panels, plot(; title="(no data)", legend=false, axis=false, border=:none))
@@ -1051,10 +1073,13 @@ function generate_plots(fort, julia, bench_dir, nn)
         j_midx = findfirst(==(m_res), j_mmodes)
         fp = (!isnothing(f_midx) && !isempty(f_bn)) ? f_bn[:, f_midx] : ComplexF64[]
         jp = (!isnothing(j_midx) && !isempty(j_bn)) ? j_bn[:, j_midx] : ComplexF64[]
-        push!(bn_panels, _profile_panel(
-            "|b_n| (norm.)", "b_n  m=$m_res (n=$nn)",
-            isempty(fp) ? Float64[] : f_psi_prof, fp, "Fortran",
-            isempty(jp) ? Float64[] : j_psi_grid, jp, "Julia"))
+        push!(
+            bn_panels,
+            _profile_panel(
+                "|b_n| (norm.)", "b_n  m=$m_res (n=$nn)",
+                isempty(fp) ? Float64[] : f_psi_prof, fp, "Fortran",
+                isempty(jp) ? Float64[] : j_psi_grid, jp, "Julia")
+        )
     end
     while length(bn_panels) < 4
         push!(bn_panels, plot(; title="(no data)", legend=false, axis=false, border=:none))
@@ -1069,10 +1094,13 @@ function generate_plots(fort, julia, bench_dir, nn)
         j_midx = findfirst(==(m_res), j_mmodes)
         fp = (!isnothing(f_midx) && !isempty(f_xigp)) ? f_xigp[:, f_midx] : ComplexF64[]
         jp = (!isnothing(j_midx) && !isempty(j_xigp)) ? j_xigp[:, j_midx] : ComplexF64[]
-        push!(xigp_panels, _profile_panel(
-            "|ξ_ψ| (norm.)", "xigradpsi  m=$m_res (n=$nn)",
-            isempty(fp) ? Float64[] : f_psi_prof, fp, "Fortran",
-            isempty(jp) ? Float64[] : j_psi_grid, jp, "Julia"))
+        push!(
+            xigp_panels,
+            _profile_panel(
+                "|ξ_ψ| (norm.)", "xigradpsi  m=$m_res (n=$nn)",
+                isempty(fp) ? Float64[] : f_psi_prof, fp, "Fortran",
+                isempty(jp) ? Float64[] : j_psi_grid, jp, "Julia")
+        )
     end
     while length(xigp_panels) < 4
         push!(xigp_panels, plot(; title="(no data)", legend=false, axis=false, border=:none))
@@ -1087,10 +1115,13 @@ function generate_plots(fort, julia, bench_dir, nn)
         j_midx = findfirst(==(m_res), j_mmodes)
         fp = (!isnothing(f_midx) && !isempty(f_xin)) ? f_xin[:, f_midx] : ComplexF64[]
         jp = (!isnothing(j_midx) && !isempty(j_xin)) ? j_xin[:, j_midx] : ComplexF64[]
-        push!(xin_panels, _profile_panel(
-            "|ξ_n| (norm.)", "xi_n  m=$m_res (n=$nn)",
-            isempty(fp) ? Float64[] : f_psi_prof, fp, "Fortran",
-            isempty(jp) ? Float64[] : j_psi_grid, jp, "Julia"))
+        push!(
+            xin_panels,
+            _profile_panel(
+                "|ξ_n| (norm.)", "xi_n  m=$m_res (n=$nn)",
+                isempty(fp) ? Float64[] : f_psi_prof, fp, "Fortran",
+                isempty(jp) ? Float64[] : j_psi_grid, jp, "Julia")
+        )
     end
     while length(xin_panels) < 4
         push!(xin_panels, plot(; title="(no data)", legend=false, axis=false, border=:none))
@@ -1105,10 +1136,13 @@ function generate_plots(fort, julia, bench_dir, nn)
         j_midx = findfirst(==(m_res), j_mmodes)
         fp = (!isnothing(f_midx) && !isempty(f_xigp1)) ? f_xigp1[:, f_midx] : ComplexF64[]
         jp = (!isnothing(j_midx) && !isempty(j_xigp1)) ? j_xigp1[:, j_midx] : ComplexF64[]
-        push!(xigp1_panels, _profile_panel(
-            "|∂ξ^ψ/∂ψ| (norm.)", "xigradpsi_dpsi  m=$m_res (n=$nn)",
-            isempty(fp) ? Float64[] : f_psi_prof, fp, "Fortran",
-            isempty(jp) ? Float64[] : j_psi_grid, jp, "Julia"))
+        push!(
+            xigp1_panels,
+            _profile_panel(
+                "|∂ξ^ψ/∂ψ| (norm.)", "xigradpsi_dpsi  m=$m_res (n=$nn)",
+                isempty(fp) ? Float64[] : f_psi_prof, fp, "Fortran",
+                isempty(jp) ? Float64[] : j_psi_grid, jp, "Julia")
+        )
     end
     while length(xigp1_panels) < 4
         push!(xigp1_panels, plot(; title="(no data)", legend=false, axis=false, border=:none))
@@ -1123,10 +1157,13 @@ function generate_plots(fort, julia, bench_dir, nn)
         j_midx = findfirst(==(m_res), j_mmodes)
         fp = (!isnothing(f_midx) && !isempty(f_xiga)) ? f_xiga[:, f_midx] : ComplexF64[]
         jp = (!isnothing(j_midx) && !isempty(j_xiga)) ? j_xiga[:, j_midx] : ComplexF64[]
-        push!(xiga_panels, _profile_panel(
-            "|ξ^α| (norm.)", "xigradalpha  m=$m_res (n=$nn)",
-            isempty(fp) ? Float64[] : f_psi_prof, fp, "Fortran",
-            isempty(jp) ? Float64[] : j_psi_grid, jp, "Julia"))
+        push!(
+            xiga_panels,
+            _profile_panel(
+                "|ξ^α| (norm.)", "xigradalpha  m=$m_res (n=$nn)",
+                isempty(fp) ? Float64[] : f_psi_prof, fp, "Fortran",
+                isempty(jp) ? Float64[] : j_psi_grid, jp, "Julia")
+        )
     end
     while length(xiga_panels) < 4
         push!(xiga_panels, plot(; title="(no data)", legend=false, axis=false, border=:none))
@@ -1157,12 +1194,12 @@ function generate_plots(fort, julia, bench_dir, nn)
     h5_path = julia["h5_path"]
     rzphi_panels = []
     for (comp_label, f_key, h5_var) in [
-            ("ξ_R",  "xi_r_fun",   "perturbed_equilibrium/response/xi_R"),
-            ("ξ_Z",  "xi_z_fun",   "perturbed_equilibrium/response/xi_Z"),
-            ("ξ_φ",  "xi_phi_fun", "perturbed_equilibrium/response/xi_phi"),
-            ("b_R",  "b_r_fun",    "perturbed_equilibrium/response/b_R"),
-            ("b_Z",  "b_z_fun",    "perturbed_equilibrium/response/b_Z"),
-            ("b_φ",  "b_phi_fun",  "perturbed_equilibrium/response/b_phi")]
+        ("ξ_R", "xi_r_fun", "perturbed_equilibrium/response/xi_R"),
+        ("ξ_Z", "xi_z_fun", "perturbed_equilibrium/response/xi_Z"),
+        ("ξ_φ", "xi_phi_fun", "perturbed_equilibrium/response/xi_phi"),
+        ("b_R", "b_r_fun", "perturbed_equilibrium/response/b_R"),
+        ("b_Z", "b_z_fun", "perturbed_equilibrium/response/b_Z"),
+        ("b_φ", "b_phi_fun", "perturbed_equilibrium/response/b_phi")]
         f_fun = get(fort, f_key, Matrix{ComplexF64}(undef, 0, 0))
 
         # Reconstruct theta-space from mode-space via modes_to_theta (applies ν phase + helicity)
@@ -1195,7 +1232,7 @@ function generate_plots(fort, julia, bench_dir, nn)
         # Real part panel
         yl_re = _percentile_ylims(f_re, j_re)
         p_re = plot(; xlabel="ψ_n", ylabel="Re($comp_label)", title="Re($comp_label) θ=0 (n=$nn)",
-                      legend=:topright, ylims=yl_re)
+            legend=:topright, ylims=yl_re)
         !isempty(f_re) && plot!(p_re, f_psi_prof, f_re; lw=2.5, color=:steelblue, label="Fortran")
         !isempty(j_re) && plot!(p_re, j_psi_grid, j_re; lw=1.0, color=:orange, label="Julia")
         push!(rzphi_panels, p_re)
@@ -1203,7 +1240,7 @@ function generate_plots(fort, julia, bench_dir, nn)
         # Imaginary part panel
         yl_im = _percentile_ylims(f_im, j_im)
         p_im = plot(; xlabel="ψ_n", ylabel="Im($comp_label)", title="Im($comp_label) θ=0 (n=$nn)",
-                      legend=:topright, ylims=yl_im)
+            legend=:topright, ylims=yl_im)
         !isempty(f_im) && plot!(p_im, f_psi_prof, f_im; lw=2.5, color=:steelblue, label="Fortran")
         !isempty(j_im) && plot!(p_im, j_psi_grid, j_im; lw=1.0, color=:orange, label="Julia")
         push!(rzphi_panels, p_im)
@@ -1228,7 +1265,7 @@ that orders-of-magnitude variation is visible.  When `rel_err=true` the colour
 scale is labelled as log₁₀(rel error).
 """
 function _mat_heatmap(M::AbstractMatrix, title_str::String, m_vals;
-                      rel_err::Bool=false, clims=nothing)
+    rel_err::Bool=false, clims=nothing)
     vals = log10.(max.(abs.(M), 1e-30))   # floor at 1e-30 to avoid -Inf
     label = rel_err ? "log₁₀(|ΔM/M|)" : "log₁₀|M|"
     kw = isnothing(clims) ? () : (clims=clims,)
@@ -1251,32 +1288,32 @@ Save `matrix_comparison_plots.png` showing element-wise |M| heatmaps for the
 fundamental response matrices L, Λ, P and the vacuum energy matrix wt0.
 
 Layout (rows × 3 columns):
-  • Col 1: |M| Fortran   Col 2: |M| Julia   Col 3: relative error |ΔM/M_fort|
-  • Row 4 (wt0): Julia only — no Fortran reference available in NetCDF outputs.
+• Col 1: |M| Fortran   Col 2: |M| Julia   Col 3: relative error |ΔM/M_fort|
+• Row 4 (wt0): Julia only — no Fortran reference available in NetCDF outputs.
 """
 function generate_matrix_plots(fort, julia, bench_dir, _nn)
     m_ctrl = fort["m_vals_ctrl"]
 
-    f_L   = fort["L_surf"]
-    f_La  = fort["Lambda"]
-    f_P   = fort["P"]
-    j_L   = julia["surface_inductance"]
-    j_La  = julia["plasma_inductance"]
-    j_P   = julia["permeability"]
+    f_L = fort["L_surf"]
+    f_La = fort["Lambda"]
+    f_P = fort["P"]
+    j_L = julia["surface_inductance"]
+    j_La = julia["plasma_inductance"]
+    j_P = julia["permeability"]
     j_wt0 = julia["wt0"]
 
     panels = Any[]
 
     for (name, f_M, j_M) in [("L  (surface inductance)", f_L, j_L),
-                               ("Λ  (plasma inductance)",  f_La, j_La),
-                               ("P  (permeability)",        f_P, j_P)]
+        ("Λ  (plasma inductance)", f_La, j_La),
+        ("P  (permeability)", f_P, j_P)]
         ok = !isempty(f_M) && !isempty(j_M) && size(f_M) == size(j_M)
         if ok
             # Use a common colour range for the first two panels
             all_vals = log10.(max.(abs.([f_M; j_M]), 1e-30))
             clims_mat = (minimum(all_vals), maximum(all_vals))
-            push!(panels, _mat_heatmap(f_M, "|$name| — Fortran",  m_ctrl; clims=clims_mat))
-            push!(panels, _mat_heatmap(j_M, "|$name| — Julia",    m_ctrl; clims=clims_mat))
+            push!(panels, _mat_heatmap(f_M, "|$name| — Fortran", m_ctrl; clims=clims_mat))
+            push!(panels, _mat_heatmap(j_M, "|$name| — Julia", m_ctrl; clims=clims_mat))
             rel_err = abs.(j_M .- f_M) ./ (abs.(f_M) .+ 1e-30)
             push!(panels, _mat_heatmap(rel_err, "Rel. error |$name|", m_ctrl; rel_err=true))
         else
@@ -1309,14 +1346,14 @@ end
 
 function main(argv=ARGS)
     opts = parse_args(argv)
-    fortran_dir      = opts.fortran_dir
-    do_plot          = opts.do_plot
-    skip_run         = opts.skip_run
+    fortran_dir = opts.fortran_dir
+    do_plot = opts.do_plot
+    skip_run = opts.skip_run
     use_file_forcing = opts.use_file_forcing
 
     isdir(fortran_dir) || error("Fortran directory not found: $fortran_dir")
     isfile(joinpath(fortran_dir, "equil.in")) || error("equil.in not found in $fortran_dir")
-    isfile(joinpath(fortran_dir, "dcon.in"))  || error("dcon.in not found in $fortran_dir")
+    isfile(joinpath(fortran_dir, "dcon.in")) || error("dcon.in not found in $fortran_dir")
 
     coil_dat_dir = joinpath(@__DIR__, "..", "src", "ForcingTerms", "coil_geometries")
 
@@ -1344,7 +1381,7 @@ function main(argv=ARGS)
     @printf "  m range: %d:%d\n" fort["mlow"] fort["mhigh"]
 
     bench_name = "fortran_vs_julia_" * basename(fortran_dir)
-    bench_dir  = joinpath(@__DIR__, bench_name)
+    bench_dir = joinpath(@__DIR__, bench_name)
     mkpath(bench_dir)
     println("\nBenchmark directory: $bench_dir")
 
@@ -1371,7 +1408,7 @@ function main(argv=ARGS)
         ratios = Float64[]
         for m in sort(m_vals)
             ja = get(julia_amps, m, NaN)
-            fa = get(fort_amps,  m, NaN)
+            fa = get(fort_amps, m, NaN)
             ratio = (!isnan(ja) && ja > 1e-40) ? fa / ja : NaN
             isnan(ratio) || push!(ratios, ratio)
             @printf "  %4d   %14.4e         %14.4e       %.4f\n" m ja fa (isnan(ratio) ? 0.0 : ratio)
@@ -1399,8 +1436,8 @@ function main(argv=ARGS)
     else
         println("  Running Julia GPEC (subprocess)...")
         isfile(h5_path) && rm(h5_path)
-        gpec_root  = joinpath(@__DIR__, "..")
-        julia_cmd  = Base.julia_cmd()
+        gpec_root = joinpath(@__DIR__, "..")
+        julia_cmd = Base.julia_cmd()
         cmd = `$julia_cmd --project=$gpec_root -e "using GeneralizedPerturbedEquilibrium; GeneralizedPerturbedEquilibrium.main([\"$bench_dir\"])"`
         run(cmd)
     end
@@ -1412,11 +1449,17 @@ function main(argv=ARGS)
 
     println()
     table_lines = build_comparison_table(fort, julia, fortran_dir, bench_dir, p.nn)
-    for line in table_lines; println(line); end
+    for line in table_lines
+        ;
+        println(line);
+    end
 
     table_path = joinpath(bench_dir, "comparison_table.txt")
     open(table_path, "w") do io
-        for line in table_lines; println(io, line); end
+        for line in table_lines
+            ;
+            println(io, line);
+        end
     end
     println("\nComparison table saved to: ", abspath(table_path))
 
diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl
index e2457c39..53e7c595 100644
--- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl
+++ b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl
@@ -129,8 +129,6 @@ shape = "nowall"
 
 [ForceFreeStates]
 local_stability_flag = true
-mat_flag = true
-ode_flag = true
 vac_flag = true
 force_termination = false
 
diff --git a/benchmarks/benchmark_diiid_kinetic_stability.jl b/benchmarks/benchmark_diiid_kinetic_stability.jl
index 0b48d0d2..9d89690c 100644
--- a/benchmarks/benchmark_diiid_kinetic_stability.jl
+++ b/benchmarks/benchmark_diiid_kinetic_stability.jl
@@ -135,8 +135,6 @@ shape = "nowall"
 
 [ForceFreeStates]
 local_stability_flag = true
-mat_flag = true
-ode_flag = true
 vac_flag = true
 force_termination = true       # Skip PE+KF post-processing — we only need FFS eigenvalues
 
diff --git a/benchmarks/gal_validation/tokamaker_beta/gal_tok_worker.jl b/benchmarks/gal_validation/tokamaker_beta/gal_tok_worker.jl
index f121f1eb..f461601e 100644
--- a/benchmarks/gal_validation/tokamaker_beta/gal_tok_worker.jl
+++ b/benchmarks/gal_validation/tokamaker_beta/gal_tok_worker.jl
@@ -24,7 +24,7 @@ function base_inputs(dir)
             "etol" => 1e-10, "force_termination" => false),
         "Wall" => Dict("shape" => "nowall", "a" => 0.2415),
         "ForceFreeStates" => Dict(
-            "local_stability_flag" => true, "mat_flag" => true, "ode_flag" => true,
+            "local_stability_flag" => true,
             "vac_flag" => true,
             "qlow" => 1.02, "qhigh" => 1e3, "sing_start" => 0,
             "nn_low" => 1, "nn_high" => 1, "delta_mlow" => 8, "delta_mhigh" => 8,
@@ -49,9 +49,13 @@ function run_pipeline(psihigh)
     eq_config = EQ.EquilibriumConfig(inputs["Equilibrium"], dir)
     equil = EQ.setup_equilibrium(eq_config, nothing)
     FFS.sing_lim!(intr, ctrl, equil)
-    xs = equil.profiles.xs; ls = zeros(length(xs), 5); FFS.compute_ballooning_stability!(ctrl, ls, equil)
+    xs = equil.profiles.xs;
+    ls = zeros(length(xs), 5);
+    FFS.compute_ballooning_stability!(ctrl, ls, equil)
     intr.locstab = cubic_interp(xs, Series(ls); extrap=ExtendExtrap())
-    intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1
+    intr.nlow = ctrl.nn_low;
+    intr.nhigh = ctrl.nn_high;
+    intr.npert = 1
     FFS.sing_find!(intr, equil)
     # Replicate main()'s post-find surface filter: drop surfaces outside [qlow, qlim].
     # Without this the worker keeps the q=1 internal-kink surface (and mis-aligns the STRIDE
@@ -81,13 +85,19 @@ end
 P = nothing
 for ph in (0.999, 0.99, 0.98, 0.97, 0.95)
     try
-        global P = run_pipeline(ph); break
+        global P = run_pipeline(ph);
+        break
     catch e
         @warn "pipeline failed at psihigh=$ph for $(basename(GEQDSK)); retrying lower" exception=e
     end
 end
 P === nothing && error("all psihigh retries failed for $(basename(GEQDSK))")
-ctrl = P.ctrl; equil = P.equil; ffit = P.ffit; intr = P.intr; vac_data = P.vac_data; dir = P.dir
+ctrl = P.ctrl;
+equil = P.equil;
+ffit = P.ffit;
+intr = P.intr;
+vac_data = P.vac_data;
+dir = P.dir
 @printf("  using psihigh=%.3f for %s\n", P.psihigh, basename(GEQDSK))
 dpm = intr.delta_prime_matrix
 m_s = [s.m[1] for s in intr.sing if s.psifac < intr.psilim]
@@ -100,7 +110,8 @@ stride_diag = [i <= ndpm ? real(dpm[i, i]) : NaN for i in 1:length(m_s)]
 open(OUTCSV, "a") do io
     for pfac in PFACS
         ctrl.gal_pfac = pfac
-        gdiag = fill(NaN, length(m_s)); gm = Int[]
+        gdiag = fill(NaN, length(m_s));
+        gm = Int[]
         try
             res = FFS.galerkin_solve(ctrl, equil, ffit, intr; vac_data=vac_data)
             gm = res.sing_m
diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml
index c55bf769..62b7930f 100644
--- a/examples/DIIID-like_gal_resistive_example/gpec.toml
+++ b/examples/DIIID-like_gal_resistive_example/gpec.toml
@@ -24,8 +24,6 @@ equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
 [ForceFreeStates]
 local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true                # Construct coefficient matrices for diagnostic purposes (required by gal_flag)
-ode_flag = true                # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                 # Edge dW scan band: dW(ψ) computed for ψ ∈ [psiedge, psilim], integration truncated at peak
@@ -55,7 +53,7 @@ dmlim               = 0.2      # Truncate integration at (last_rational_q + dmli
 
 # Outer-region singular Galerkin Δ′ solver with rpec coil columns.
 # This case exercises the gal solve AND the rpec coil/edge path (delta_coil) for regression tracking.
-gal_flag               = true  # Enable the outer-region singular Galerkin Δ′ solve (requires mat_flag = true)
+gal_flag               = true  # Enable the outer-region singular Galerkin Δ′ solve
 gal_solver             = "LU"  # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU".
 gal_nx                 = 256   # Hermite-cubic elements per interval between singular surfaces
 gal_nq                 = 6     # Gauss-Lobatto quadrature order per element
diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
index ffc03cb5..a6c6612a 100644
--- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
+++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
@@ -24,8 +24,6 @@ equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
 [ForceFreeStates]
 local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true                # Construct coefficient matrices for diagnostic purposes (required by gal_flag)
-ode_flag = true                # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                 # Edge dW scan band: dW(ψ) computed for ψ ∈ [psiedge, psilim], integration truncated at peak
@@ -55,7 +53,7 @@ dmlim               = 0.2      # Truncate integration at (last_rational_q + dmli
 
 # Outer-region singular Galerkin Δ′ solver with rpec coil columns.
 # This case exercises the gal solve AND the rpec coil/edge path (delta_coil) for regression tracking.
-gal_flag               = true  # Enable the outer-region singular Galerkin Δ′ solve (requires mat_flag = true)
+gal_flag               = true  # Enable the outer-region singular Galerkin Δ′ solve
 gal_solver             = "LU"  # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU".
 gal_nx                 = 256   # Hermite-cubic elements per interval between singular surfaces
 gal_nq                 = 6     # Gauss-Lobatto quadrature order per element
diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml
index 5c7be179..a36d2127 100644
--- a/examples/DIIID-like_ideal_example/gpec.toml
+++ b/examples/DIIID-like_ideal_example/gpec.toml
@@ -29,8 +29,6 @@ equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
 [ForceFreeStates]
 local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true                # Construct coefficient matrices for diagnostic purposes
-ode_flag = true                # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
diff --git a/examples/DIIID-like_ideal_example_IMAS/gpec.toml b/examples/DIIID-like_ideal_example_IMAS/gpec.toml
index fec92d8b..0893460a 100644
--- a/examples/DIIID-like_ideal_example_IMAS/gpec.toml
+++ b/examples/DIIID-like_ideal_example_IMAS/gpec.toml
@@ -24,8 +24,6 @@ shape = "nowall"               # Wall shape (nowall, conformal, elliptical, dee,
 nn_low               = 1       # Smallest toroidal mode number to include
 nn_high              = 1       # Largest toroidal mode number to include
 vac_flag             = true    # Compute plasma, vacuum, and total energies for free-boundary modes
-ode_flag             = true    # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
-mat_flag             = true    # Construct coefficient matrices for diagnostic purposes
 local_stability_flag = false   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 write_outputs_to_HDF5 = false  # Write outputs to the HDF5 file
 verbose              = false   # Enable verbose logging
diff --git a/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl b/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl
index 178da6ad..df964153 100644
--- a/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl
+++ b/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl
@@ -26,41 +26,42 @@ const GEQDSK = joinpath(@__DIR__, "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmo
 # Only the [Equilibrium] section differs between IMAS and g-file paths.
 const SHARED_CONFIG = Dict(
     "ForceFreeStates" => Dict(
-        "nn_low"               => 1,
-        "nn_high"              => 1,
-        "vac_flag"             => true,
-        "ode_flag"             => true,
-        "mat_flag"             => true,
+        "nn_low" => 1,
+        "nn_high" => 1,
+        "vac_flag" => true,
         "local_stability_flag" => false,
-        "write_outputs_to_HDF5"=> false,
-        "verbose"              => false,
-        "delta_mhigh"          => 16,
-        "qlow"                 => 1.02,
+        "write_outputs_to_HDF5" => false,
+        "verbose" => false,
+        "delta_mhigh" => 16,
+        "qlow" => 1.02
     ),
-    "Wall" => Dict("shape" => "nowall"),
+    "Wall" => Dict("shape" => "nowall")
 )
 
 const EQ_COMMON = Dict(
-    "jac_type"   => "boozer",
-    "psilow"     => 0.01,
-    "psihigh"    => 0.9995,
-    "mpsi"       => 128,
-    "mtheta"     => 256,
+    "jac_type" => "boozer",
+    "psilow" => 0.01,
+    "psihigh" => 0.9995,
+    "mpsi" => 128,
+    "mtheta" => 256
 )
 
 @info "Loading g-file and converting to IMAS (COCOS 11)..."
-g  = EFIT.readg(GEQDSK)
+g = EFIT.readg(GEQDSK)
 dd = IMASdd.dd()
 EFIT.geqdsk2imas!([g], dd; dd_cocos=11)
 
 @info "Running GPEC from IMAS dd..."
-tmpdir_imas = mktempdir(prefix="gpec_imas_")
+tmpdir_imas = mktempdir(; prefix="gpec_imas_")
 try
     config_imas = merge(SHARED_CONFIG, Dict("Equilibrium" =>
         merge(EQ_COMMON, Dict("eq_type" => "imas", "eq_filename" => "from_dd", "imas_cocos" => 11))))
-    open(joinpath(tmpdir_imas, "gpec.toml"), "w") do io; TOML.print(io, config_imas); end
+    open(joinpath(tmpdir_imas, "gpec.toml"), "w") do io
+        ;
+        TOML.print(io, config_imas);
+    end
     result_imas = GPEC.main([tmpdir_imas]; dd=dd)
-    global et_imas  = real(result_imas.vac_data.et[1])
+    global et_imas = real(result_imas.vac_data.et[1])
     global mpert_imas = result_imas.intr.mpert
     GPEC.write_imas(dd, result_imas)
     @assert dd.mhd_linear.time_slice[1].toroidal_mode[1].energy_perturbed ≈ et_imas
@@ -70,13 +71,16 @@ end
 @info @sprintf "IMAS path:   δW(n=1) = %.6f  [mpert = %d]" et_imas mpert_imas
 
 @info "Running GPEC from g-file (reference)..."
-tmpdir_gfile = mktempdir(prefix="gpec_gfile_")
+tmpdir_gfile = mktempdir(; prefix="gpec_gfile_")
 try
     config_gfile = merge(SHARED_CONFIG, Dict("Equilibrium" =>
         merge(EQ_COMMON, Dict("eq_type" => "efit", "eq_filename" => GEQDSK))))
-    open(joinpath(tmpdir_gfile, "gpec.toml"), "w") do io; TOML.print(io, config_gfile); end
+    open(joinpath(tmpdir_gfile, "gpec.toml"), "w") do io
+        ;
+        TOML.print(io, config_gfile);
+    end
     result_gfile = GPEC.main([tmpdir_gfile])
-    global et_gfile  = real(result_gfile.vac_data.et[1])
+    global et_gfile = real(result_gfile.vac_data.et[1])
     global mpert_gfile = result_gfile.intr.mpert
 finally
     rm(tmpdir_gfile; recursive=true)
@@ -89,7 +93,7 @@ println("\n" * "="^60)
 println("  IMAS vs g-file comparison (n=1, nowall, boozer)")
 println("="^60)
 @printf "  G-file δW(n=1) : %+.6f  [mpert = %d]\n" et_gfile mpert_gfile
-@printf "  IMAS   δW(n=1) : %+.6f  [mpert = %d]\n" et_imas  mpert_imas
+@printf "  IMAS   δW(n=1) : %+.6f  [mpert = %d]\n" et_imas mpert_imas
 @printf "  Relative diff  : %.2f%%\n" rel_diff * 100
 println("="^60)
 
diff --git a/examples/LAR_beta_scan/gpec.toml b/examples/LAR_beta_scan/gpec.toml
index e2f7a924..6e5c664f 100644
--- a/examples/LAR_beta_scan/gpec.toml
+++ b/examples/LAR_beta_scan/gpec.toml
@@ -35,8 +35,6 @@ a     = 20                     # Distance from plasma (conformal) or shape param
 
 [ForceFreeStates]
 local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true                # Construct coefficient matrices for diagnostic purposes
-ode_flag = true                # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 qlow         = 1.02            # Integration initiated at q determined by min(q0, qlow)
diff --git a/examples/LAR_epsilon_scan/gpec.toml b/examples/LAR_epsilon_scan/gpec.toml
index c1b195e4..b5cc871e 100644
--- a/examples/LAR_epsilon_scan/gpec.toml
+++ b/examples/LAR_epsilon_scan/gpec.toml
@@ -36,8 +36,6 @@ a     = 20                     # Distance from plasma (conformal) or shape param
 
 [ForceFreeStates]
 local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true                # Construct coefficient matrices for diagnostic purposes
-ode_flag = true                # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 qlow         = 1.02            # Integration initiated at q determined by min(q0, qlow)
diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml
index 77b48f29..82367adc 100644
--- a/examples/Solovev_ideal_example/gpec.toml
+++ b/examples/Solovev_ideal_example/gpec.toml
@@ -46,8 +46,6 @@ reg_spot = 0.05                         # Regularization width for singular surf
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml
index 017428d6..991cfc4d 100644
--- a/examples/Solovev_ideal_example_3D/gpec.toml
+++ b/examples/Solovev_ideal_example_3D/gpec.toml
@@ -16,8 +16,6 @@ force_termination = false               # Terminate after equilibrium setup (ski
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 1.0                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
diff --git a/examples/Solovev_ideal_example_multi_n/gpec.toml b/examples/Solovev_ideal_example_multi_n/gpec.toml
index 09a24132..3fb3f7b0 100644
--- a/examples/Solovev_ideal_example_multi_n/gpec.toml
+++ b/examples/Solovev_ideal_example_multi_n/gpec.toml
@@ -28,8 +28,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
diff --git a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml
index 8853b817..8b2b7500 100644
--- a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml
+++ b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml
@@ -27,8 +27,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
diff --git a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml
index 677963f2..c7a51c72 100644
--- a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml
+++ b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml
@@ -27,8 +27,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml
index ecc36a6a..f1d3ee96 100644
--- a/examples/Solovev_kinetic_NTV_example/gpec.toml
+++ b/examples/Solovev_kinetic_NTV_example/gpec.toml
@@ -48,8 +48,6 @@ reg_spot = 0.05                         # Regularization width for singular surf
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml
index 15ce5854..ac27543c 100644
--- a/examples/Solovev_kinetic_calculated_example/gpec.toml
+++ b/examples/Solovev_kinetic_calculated_example/gpec.toml
@@ -27,8 +27,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = false              # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW scan band: dW(ψ) computed for ψ ∈ [psiedge, psilim], integration truncated at peak
diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml
index 66ace366..d3f61b39 100644
--- a/examples/a10_kinetic_example/gpec.toml
+++ b/examples/a10_kinetic_example/gpec.toml
@@ -23,8 +23,6 @@ shape = "nowall"               # Wall shape (nowall, conformal, elliptical, dee,
 
 [ForceFreeStates]
 local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true                # Construct coefficient matrices for diagnostic purposes
-ode_flag = true                # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 set_psilim_via_dmlim = false   # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge
diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl
index 249e6bd7..1bbee56b 100644
--- a/src/Analysis/ForceFreeStates.jl
+++ b/src/Analysis/ForceFreeStates.jl
@@ -41,7 +41,7 @@ function plot_mode_displacement(h5path; modes=1:5, save_path=nothing)
     mhigh = mlow + mpert - 1
     dW = isempty(et) ? nothing : et[1]
     title_str = isnothing(dW) ? "Least stable mode" :
-        "Least stable mode, δW = $(round(real(dW), sigdigits=4))"
+                "Least stable mode, δW = $(round(real(dW), sigdigits=4))"
 
     p = plot(;
         xlims=(0, 1),
@@ -425,7 +425,7 @@ end
 
 Plot `cond(F̄)` vs ψ from the kinetic-singular-surface scan stored in
 `singular/kinetic/` (populated when ForceFreeStates runs with
-`kinetic_factor > 0`, `ode_flag = true`, `singfac_min > 0`).
+`kinetic_factor > 0`, `singfac_min > 0`).
 
 `F̄` is the kinetic Euler-Lagrange matrix formed by Schur-reducing the six
 kinetic matrices against the ideal A/B/C/D/E/H blocks (Logan 2015 Appendix
@@ -471,7 +471,7 @@ function plot_cond_fbar(h5path; save_path=nothing, zoom=false)
     end
 
     if isempty(scan_psi)
-        @warn "No kinetic-singular-surface scan in $h5path — rerun with kinetic_factor>0, ode_flag=true, singfac_min>0"
+        @warn "No kinetic-singular-surface scan in $h5path — rerun with kinetic_factor>0, singfac_min>0"
         return nothing
     end
 
diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl
index 025007ee..7f727ea5 100644
--- a/src/ForceFreeStates/ForceFreeStatesStructs.jl
+++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl
@@ -205,8 +205,6 @@ A mutable struct containing control parameters for stability analysis, set by th
 
   - `verbose::Bool` - Enable verbose output
   - `local_stability_flag::Bool` - Enable local stability analysis (`D_I` and ballooning)
-  - `mat_flag::Bool` - Enable matrix output
-  - `ode_flag::Bool` - Enable ODE integration diagnostics
   - `vac_flag::Bool` - Enable vacuum region calculation
   - `mthvac::Int` - Number of vacuum poloidal grid points (corresponds to `mtheta` in VacuumInput)
   - `nzvac::Int` - Number of vacuum toroidal grid points (corresponds to `nzeta` in VacuumInput3D)
@@ -222,7 +220,6 @@ A mutable struct containing control parameters for stability analysis, set by th
   - `numsteps_init::Int` - Initial array size for ODE data storage
   - `numunorms_init::Int` - Initial array size for solution normalization data
   - `singfac_min::Float64` - Fractional distance from rational q at which ideal jump condition is enforced
-  - `cyl_flag::Bool` - Make delta_mlow and delta_mhigh set the actual m truncation bounds. Default is to expand (n*qmin-4, n*qmax).
   - `set_psilim_via_dmlim::Bool` - Truncate the integration domain at `(last_rational_q + dmlim) / n` rather than at `qhigh` / `psihigh`. Fortran STRIDE found that truncating ~20 % above the outermost rational (`dmlim = 0.2`) avoids a numerical kink instability in δW that appears when the integration ends too close to or just below a rational surface. **For diverted equilibria where q → ∞ at the separatrix** (e.g. DIII-D geqdsks, the bulk of production use) this costs negligible physical domain because rationals get arbitrarily dense near the LCFS — `set_psilim_via_dmlim = true` is the safe and recommended default. **For limited circular / analytical equilibria with finite q at the edge** (Solovev, LAR scans), rationals are sparse and 20 % above the last rational chops off too much edge, so set `set_psilim_via_dmlim = false` and let `qhigh` / `psihigh` control the truncation. Multi-`n` runs are not supported by this truncation (the "outermost rational + dmlim / n" depends on which `n`); when `set_psilim_via_dmlim = true` with `nn_low != nn_high`, `sing_lim!` warns and falls back to `qhigh` / `psihigh`. Default `true`.
   - `dmlim::Float64` - Distance beyond last rational surface (normalised ∈ [0,1) in units of 1/n). Only used when `set_psilim_via_dmlim` is true. Fortran STRIDE convention is 0.2 (truncate 20 % of one rational-surface spacing above the last surface), retained here.
   - `sing_order::Int` - Order of singular layer (Frobenius) expansion at rational surfaces. Default 6 (Fortran STRIDE convention for Δ' calculations; lower values trade accuracy for speed).
@@ -248,8 +245,6 @@ A mutable struct containing control parameters for stability analysis, set by th
 @kwdef mutable struct ForceFreeStatesControl
     verbose::Bool = true
     local_stability_flag::Bool = false
-    mat_flag::Bool = false
-    ode_flag::Bool = false
     vac_flag::Bool = false
     mthvac::Int = 480
     nzvac::Int = 1
@@ -265,7 +260,6 @@ A mutable struct containing control parameters for stability analysis, set by th
     numsteps_init::Int = 4000
     numunorms_init::Int = 100
     singfac_min::Float64 = 1e-4   # Matches Fortran STRIDE; required nonzero for use_parallel path.
-    cyl_flag::Bool = false
     set_psilim_via_dmlim::Bool = true   # Safe default for diverted equilibria (most production use); set false for limited/analytical (LAR, Solovev). Auto-skipped for multi-n. See docstring.
     dmlim::Float64 = 0.2
     sing_order::Int = 6
@@ -342,7 +336,7 @@ end
     kmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts)
     gmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts)
 
-    # Ideal A,B,C splines preserved before kinetic overwrite (for mat_flag output)
+    # Ideal A,B,C splines preserved before kinetic overwrite
     amats_ideal::S = _empty_series_interp_complex(numpert_total^2, itp_opts)
     bmats_ideal::S = _empty_series_interp_complex(numpert_total^2, itp_opts)
     cmats_ideal::S = _empty_series_interp_complex(numpert_total^2, itp_opts)
diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl
index 091855e5..52acb9ba 100644
--- a/src/ForceFreeStates/Free.jl
+++ b/src/ForceFreeStates/Free.jl
@@ -49,9 +49,7 @@ and data dumping.
     dV_dpsi = equil.profiles.dVdpsi_spline(psilim)
 
     # Compute plasma response matrix W = U₂ * U₁⁻¹
-    if ctrl.ode_flag
-        @views wp .= (odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2
-    end
+    @views wp .= (odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2
 
     # Compute vacuum response matrix in-place (handles 2D single-n, 2D multi-n block-diagonal, and 3D)
     vac_inputs = Vacuum.VacuumInput(equil, psilim, ctrl.mthvac, ctrl.nzvac, mlow:mhigh, nlow:nhigh)
diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl
index 389c6771..a1434b85 100644
--- a/src/ForceFreeStates/Kinetic.jl
+++ b/src/ForceFreeStates/Kinetic.jl
@@ -26,7 +26,7 @@ function make_kinetic_matrix(
     ffit::FourFitVars,
     intr::ForceFreeStatesInternal,
     metric::MetricData;
-    calculated_source::Union{Nothing,Function}=nothing,
+    calculated_source::Union{Nothing,Function}=nothing
 )
     xs = metric.xs
     mpsi = length(xs)
@@ -248,7 +248,7 @@ function _compute_fkg_matrices!(
     ffit.r3mats = cubic_interp(xs, Series(r3_flat); ffit.itp_opts...)
     ffit.gaats = cubic_interp(xs, Series(ga_flat); ffit.itp_opts...)
 
-    # Preserve ideal A/B/C splines before overwrite (for mat_flag output)
+    # Preserve ideal A/B/C splines before overwrite
     ffit.amats_ideal = ffit.amats
     ffit.bmats_ideal = ffit.bmats
     ffit.cmats_ideal = ffit.cmats
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index edee8303..2e10a1bd 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -59,7 +59,7 @@ using .ForceFreeStates: find_kinetic_singular_surfaces!
 using .ForceFreeStates: eulerlagrange_integration, free_run!
 using .ForceFreeStates: galerkin_solve, write_galerkin!, GalerkinResult, gal_matched_odestate
 
-const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry")
+const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag")
 const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc")
 
 # Drop deprecated keys from a parsed gpec.toml section so legacy files keep parsing
@@ -370,10 +370,7 @@ function main_from_inputs(
     if ctrl.delta_mlow < 0 || ctrl.delta_mhigh < 0
         error("Negative delta_mlow or delta_mhigh not allowed")
     end
-    if ctrl.cyl_flag
-        intr.mlow = ctrl.delta_mlow
-        intr.mhigh = ctrl.delta_mhigh
-    elseif ctrl.sing_start == 0
+    if ctrl.sing_start == 0
         intr.mlow = trunc(Int, min(intr.nlow * equil.params.qmin, 0)) - 4 - ctrl.delta_mlow
         intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
     else
@@ -388,64 +385,54 @@ function main_from_inputs(
     intr.numpert_total = intr.mpert * intr.npert
 
     # Fit equilibrium quantities to Fourier-spline functions.
-    if ctrl.mat_flag || ctrl.ode_flag
-        if ctrl.verbose
-            @info "Run parameters:\n" *
-                  "   q0 = $(@sprintf("%.3f", equil.params.q0)), qmin = $(@sprintf("%.3f", equil.params.qmin)), qmax = $(@sprintf("%.3f", equil.params.qmax)), q95 = $(@sprintf("%.3f", equil.params.q95))\n" *
-                  "   qlim = $(@sprintf("%.3f", intr.qlim)), psilim = $(@sprintf("%.3f", intr.psilim))\n" *
-                  "   betat = $(@sprintf("%.3f", equil.params.betat)), betan = $(@sprintf("%.3f", equil.params.betan)), betap1 = $(@sprintf("%.3f", equil.params.betap1))\n" *
-                  "   mlow = $(@sprintf("%4i", intr.mlow)), mhigh = $(@sprintf("%4i", intr.mhigh)), mpert = $(@sprintf("%4i", intr.mpert))\n" *
-                  "   nlow = $(@sprintf("%4i", intr.nlow)), nhigh = $(@sprintf("%4i", intr.nhigh)), npert = $(@sprintf("%4i", intr.npert))"
-        end
+    if ctrl.verbose
+        @info "Run parameters:\n" *
+              "   q0 = $(@sprintf("%.3f", equil.params.q0)), qmin = $(@sprintf("%.3f", equil.params.qmin)), qmax = $(@sprintf("%.3f", equil.params.qmax)), q95 = $(@sprintf("%.3f", equil.params.q95))\n" *
+              "   qlim = $(@sprintf("%.3f", intr.qlim)), psilim = $(@sprintf("%.3f", intr.psilim))\n" *
+              "   betat = $(@sprintf("%.3f", equil.params.betat)), betan = $(@sprintf("%.3f", equil.params.betan)), betap1 = $(@sprintf("%.3f", equil.params.betap1))\n" *
+              "   mlow = $(@sprintf("%4i", intr.mlow)), mhigh = $(@sprintf("%4i", intr.mhigh)), mpert = $(@sprintf("%4i", intr.mpert))\n" *
+              "   nlow = $(@sprintf("%4i", intr.nlow)), nhigh = $(@sprintf("%4i", intr.nhigh)), npert = $(@sprintf("%4i", intr.npert))"
+    end
 
-        # Compute metric tensor
-        metric = make_metric(equil, intr.mpert)
+    # Compute metric tensor
+    metric = make_metric(equil, intr.mpert)
 
-        if ctrl.verbose
-            @info "Computing F, G, and K matrices"
-        end
+    if ctrl.verbose
+        @info "Computing F, G, and K matrices"
+    end
 
-        # Compute matrices and populate FourFitVars struct
-        ffit = make_matrix(equil, intr, metric)
+    # Compute matrices and populate FourFitVars struct
+    ffit = make_matrix(equil, intr, metric)
 
-        if ctrl.kinetic_factor > 0
-            if ctrl.verbose
-                @info "Computing kinetic matrices (source: $(ctrl.kinetic_source), factor: $(ctrl.kinetic_factor))"
-            end
-            # Inject the KineticForces callback so the "calculated" source can
-            # invoke compute_calculated_kinetic_matrices without ForceFreeStates
-            # importing KineticForces (which would invert the load order).
-            calculated_cb = (c, e, i, m, f) ->
-                KineticForces.compute_calculated_kinetic_matrices(
-                    c, e, i, m, f;
-                    kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles)
-            make_kinetic_matrix(ctrl, equil, ffit, intr, metric;
-                calculated_source=calculated_cb)
-
-            # Find kinetically-displaced singular surfaces (zeros of det(F̄)) for ODE crossings.
-            # Matches Fortran ksing_find (sing.f:1486-1616). singfac_min > 0 gates crossings;
-            # singfac_min == 0 preserves single-chunk behavior.
-            if ctrl.ode_flag && ctrl.singfac_min > 0
-                find_kinetic_singular_surfaces!(ffit, equil, intr)
-            end
+    if ctrl.kinetic_factor > 0
+        if ctrl.verbose
+            @info "Computing kinetic matrices (source: $(ctrl.kinetic_source), factor: $(ctrl.kinetic_factor))"
+        end
+        # Inject the KineticForces callback so the "calculated" source can
+        # invoke compute_calculated_kinetic_matrices without ForceFreeStates
+        # importing KineticForces (which would invert the load order).
+        calculated_cb = (c, e, i, m, f) ->
+            KineticForces.compute_calculated_kinetic_matrices(
+                c, e, i, m, f;
+                kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles)
+        make_kinetic_matrix(ctrl, equil, ffit, intr, metric;
+            calculated_source=calculated_cb)
+
+        # Find kinetically-displaced singular surfaces (zeros of det(F̄)) for ODE crossings.
+        # Matches Fortran ksing_find (sing.f:1486-1616). singfac_min > 0 gates crossings;
+        # singfac_min == 0 preserves single-chunk behavior.
+        if ctrl.singfac_min > 0
+            find_kinetic_singular_surfaces!(ffit, equil, intr)
         end
-
-        # NOTE: Asymptotic calculations for ideal ForceFreeStates are now computed on-demand during
-        # singular surface crossings in cross_ideal_singular_surf!. This makes it clear that
-        # asymptotics are only needed for ideal ForceFreeStates and are not inherent properties of
-        # the singular surface.
-
     end
 
     # Integrate Euler-Lagrange Equation
-    if ctrl.ode_flag
-        if ctrl.verbose
-            @info "Integrating Euler-Lagrange equation"
-        end
-        odet, fm_propagators, fm_chunks, fm_S_left = eulerlagrange_integration(ctrl, equil, ffit, intr)
-        if odet.nzero > 0 && ctrl.verbose
-            @warn "Fixed-boundary mode unstable for n = $nstring"
-        end
+    if ctrl.verbose
+        @info "Integrating Euler-Lagrange equation"
+    end
+    odet, fm_propagators, fm_chunks, fm_S_left = eulerlagrange_integration(ctrl, equil, ffit, intr)
+    if odet.nzero > 0 && ctrl.verbose
+        @warn "Fixed-boundary mode unstable for n = $nstring"
     end
 
     # Compute free boundary energies
@@ -482,7 +469,6 @@ function main_from_inputs(
     gal_data = nothing
     if ctrl.gal_flag
         gal_start = time()
-        ctrl.mat_flag || error("gal_flag=true requires mat_flag=true (needs the F/G/K matrix splines)")
         gal_data = galerkin_solve(ctrl, equil, ffit, intr; vac_data=ctrl.vac_flag ? vac_data : nothing)
         @info "Galerkin solve completed in $(@sprintf("%.3f", time() - gal_start)) s"
     end
@@ -745,7 +731,6 @@ function write_outputs_to_HDF5(
         out_h5["locstab/alpha_critical"] = ballooning_boundary.alpha_critical
 
         # Write integration data
-        # TODO: technically this should only be written if ode_flag is true, but that's going to get deprecated eventually
         out_h5["integration/nstep"] = odet.step            # Number of saved solution snapshots
         out_h5["integration/nstep_total"] = odet.total_steps  # Total ODE solver steps taken
         out_h5["integration/psi"] = odet.psi_store
@@ -840,8 +825,8 @@ function write_outputs_to_HDF5(
             out_h5["kinetic/kinetic_factor"] = ctrl.kinetic_factor
         end
 
-        # Write fundamental matrices on the ψ grid when mat_flag is enabled
-        if ctrl.mat_flag && ffit !== nothing
+        # Write fundamental matrices on the ψ grid
+        if ffit !== nothing
             xs = equil.rzphi_xs
             npsi = length(xs)
             np = intr.numpert_total
diff --git a/test/test_data/regression_solovev_ideal_example/gpec.toml b/test/test_data/regression_solovev_ideal_example/gpec.toml
index 6608c382..e6a97d2e 100644
--- a/test/test_data/regression_solovev_ideal_example/gpec.toml
+++ b/test/test_data/regression_solovev_ideal_example/gpec.toml
@@ -26,8 +26,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
diff --git a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml
index 90c8cf56..0f0bc5c4 100644
--- a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml
+++ b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml
@@ -26,8 +26,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
diff --git a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml
index c6fe1ab2..5b87267a 100644
--- a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml
+++ b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml
@@ -27,8 +27,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = false              # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
diff --git a/test/test_data/regression_solovev_kinetic_example/gpec.toml b/test/test_data/regression_solovev_kinetic_example/gpec.toml
index 763ddf41..4d8f0310 100644
--- a/test/test_data/regression_solovev_kinetic_example/gpec.toml
+++ b/test/test_data/regression_solovev_kinetic_example/gpec.toml
@@ -26,8 +26,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = false              # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
diff --git a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml
index 189f13a2..b11caab8 100644
--- a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml
+++ b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml
@@ -26,8 +26,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = false              # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
diff --git a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml
index 6d3b5b9b..3aab7516 100644
--- a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml
+++ b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml
@@ -22,8 +22,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = false              # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW scan band: dW(ψ) computed for ψ ∈ [psiedge, psilim], integration truncated at peak