function F = energy_env_steady_equations_lsqnonlin(x)

    % Unpack variables
    k_fossil  = x(1);
    k_rest    = x(2);
    k_total   = x(3);
    l         = x(4);
    inv_fossil = x(5);
    inv_rest   = x(6);
    cap       = x(7);
    c         = x(8);
    r         = x(9);
    w         = x(10);
    a         = x(11);
    z         = x(12);
    tau       = x(13);
    T_tax     = x(14);
    E         = x(15);
    E_fossil  = x(16);
    E_renew   = x(17);
    F         = x(18);
    ND        = x(19);
    X         = x(20);
    Z         = x(21);
    y         = x(22);
    g         = x(23);
    G         = x(24);
    A_x       = x(25);
    d         = x(26);

    % Parameters
    beta = 0.99;
    alpha = 0.33;
    delta = 0.025;
    gamma = 0.05;
    eta = 2;
    phi = 0.1;
    chi = 0.2;
    psi = 0.15;
    lambda = 0.1;
    a1 = 0.3;
    a2 = 0.6;
    phi_e = 2;
    mu = 0.5;
    epsilon_E = 0.8;
    a_fossil = 0.4;
    a_renew = 0.4;
    A_F = 1;
    A_R = 1;
    delta_Z = 0.05;
    delta_F = 0.03;
    n = 0.5;
    theta = 0.7;

    % Model equations
    F = zeros(26, 1);

    F(1) = 1/c - beta * (1/c) * (1 + r);
    F(2) = l^phi_e - w / c;
    F(3) = c + inv_fossil + inv_rest + T_tax - w*l - r*(inv_fossil + inv_rest) - d;
    F(4) = k_fossil - (1 - delta) * k_fossil - inv_fossil;
    F(5) = k_rest - (1 - delta) * k_rest - inv_rest;
    F(6) = k_total - (k_fossil + k_rest);
    F(7) = E - (mu * E_renew^((epsilon_E - 1)/epsilon_E) + (1 - mu) * E_fossil^((epsilon_E - 1)/epsilon_E))^(epsilon_E / (epsilon_E - 1));
    F(8) = E_fossil - A_F * k_fossil^a_fossil * l^(1 - a_fossil) * F;
    F(9) = E_renew - A_R * k_rest^a_renew * l^(1 - a_renew);
    F(10)= X - A_x * n * E;
    F(11)= Z - X / delta_Z;
    F(12)= F - ND / delta_F;
    F(13)= y - exp(a) * exp(-gamma * Z) * k_total^a1 * l^a2 * E^(1 - a1 - a2);
    F(14)= g - eta * exp(z) * exp(a) * k_total^alpha * (1 - lambda * cap);
    F(15)= G - theta * g;
    F(16)= r - a1 * y / k_total;
    F(17)= w - a2 * y / l;
    F(18)= d;
    F(19)= y - c - E - inv_fossil - inv_rest - G;
    F(20)= a;
    F(21)= z;
    F(22)= tau - 0.02;
    F(23)= T_tax - 0.015;
    F(24)= cap - 0.02;
    F(25)= A_x - 0.1;
    F(26)= ND - 0.03;

  % Generate Dynare-ready initval block
    var_names = {
        'c','l','inv_fossil','inv_rest','T_tax','w','r','d','a','z','tau','cap','g','Z',
        'F','X','E_fossil','E_renew','y','k_fossil','k_rest','G','A_x','ND','k_total','E'};

    fprintf('\n\nPaste this into your Dynare .mod file:\n');
    fprintf('initval;\n');
    for i = 1:length(var_names)
        fprintf('%-12s = %.16f;\n', var_names{i}, x(i));
    end
    fprintf('end;\n');
end
