
function [ss_values, resnorm, residuals, exitflag] = run_lsqnonlin_energy_env_L()

    % Initial guess based on steady_state_model values (in levels)
    x0 = [
        1.5;  % k_fossil
        1.5;  % k_rest
        3.0;  % k_total
        0.7;  % l
        0.25; % inv_fossil
        0.25; % inv_rest
        0.1;  % cap
        4.5;  % c
        0.05; % tau
        0.4;  % T_tax
        0;    % a
        0;    % z
        0.1;  % A_x
        0.1;  % ND
        1.5;  % F
        1.2;  % E
        0.3;  % X
        8.5;  % Z
        0.3;  % E_fossil
        1.2;  % E_renew
        0.8;  % y
        8.5;  % g
        1.05; % G
        1.0;  % r
        0.2;  % w
        1.5;  % d
    ];

    % Options
    options = optimoptions('lsqnonlin', 'Display', 'iter-detailed', 'FunctionTolerance', 1e-8, 'StepTolerance', 1e-8, 'OptimalityTolerance', 1e-8, 'MaxIterations', 3000, 'MaxFunctionEvaluations', 8000);

    % Run solver
    [ss_values, resnorm, residuals, exitflag] = lsqnonlin(@energy_env_steady_equations, x0, [], [], options);

    % Display results
    disp("✅ Steady State Values:");
    disp(ss_values');
    disp("📉 Function Residuals:");
    disp(residuals);
    disp("🚦 Exit Flag:");
    disp(exitflag);
end
