function [ys, params, check] = example_PissaridesII_steadystate(ys, exo, M_, options_)
    % Initialize check
    check = 0;

    % Retrieve parameters by direct assignment
    c     = M_.params(strmatch('c',     M_.param_names, 'exact'));
    delta = M_.params(strmatch('delta', M_.param_names, 'exact'));
    beta  = M_.params(strmatch('beta',  M_.param_names, 'exact'));
    mu    = M_.params(strmatch('mu',    M_.param_names, 'exact'));
    alpha = M_.params(strmatch('alpha', M_.param_names, 'exact'));
    eta   = M_.params(strmatch('eta',   M_.param_names, 'exact'));
    b     = M_.params(strmatch('b',     M_.param_names, 'exact'));
    rho   = M_.params(strmatch('rho',   M_.param_names, 'exact'));

    % Retrieve exogenous variables
    epsilon = 0; % Default value if exo is empty
    if length(exo) >= 1
        epsilon = exo(1);
    end

    % Compute steady-state productivity Z
    Z = exp(epsilon / (1 - rho));

    % Compute A
    A = 1 - beta * (1 - delta);

    % Define function to solve for theta and n simultaneously
    % Set up a system of two equations:
    % 1. Free entry condition
    % 2. Employment steady-state condition

    % Initial guesses for theta and n
    x0 = [1; 0.9];

    % Solve the system numerically
    options_fsolve = optimset('Display', 'off');
    [x_sol, fval, exitflag] = fsolve(@steady_state_system, x0, options_fsolve);

    if exitflag <= 0
        check = 1;
        return;
    end

    theta = x_sol(1);
    n = x_sol(2);

    % Wage
    w = eta * (Z + c * theta) + (1 - eta) * b;

    % Job Value
    J = (Z - w) / A;

    % Value of being employed and unemployed
    f = mu * theta^(1 - alpha);
    A_w = 1 - beta * (1 - delta);
    A_u = 1 - beta * (1 - f);

    % Set up the system of equations to solve for W and U
    % Bellman equations at steady state:
    % A_w * W - beta * delta * U = w
    % A_u * U - beta * f * W = b

    % Coefficient matrix
    coeff_matrix = [A_w, -beta * delta;
                    -beta * f,   A_u];

    % Right-hand side vector
    rhs = [w; b];

    % Solve for W and U
    sol = coeff_matrix \ rhs;
    W = sol(1);
    U = sol(2);

    % Since V = 0 at steady state
    V = 0;

    % Compute total output Y
    Y = Z * n;

    % Prepare outputs
    ys = zeros(M_.endo_nbr, 1);

    % Assign steady-state values
    ys(strcmp(M_.endo_names, 'Z'))     = Z;
    ys(strcmp(M_.endo_names, 'theta')) = theta;
    ys(strcmp(M_.endo_names, 'n'))     = n;
    ys(strcmp(M_.endo_names, 'w'))     = w;
    ys(strcmp(M_.endo_names, 'J'))     = J;
    ys(strcmp(M_.endo_names, 'V'))     = V;
    ys(strcmp(M_.endo_names, 'W'))     = W;
    ys(strcmp(M_.endo_names, 'U'))     = U;
    ys(strcmp(M_.endo_names, 'Y'))     = Y;

    % Assign params (optional)
    params = M_.params;

    % Nested function defining the system of equations
    function F = steady_state_system(x)
        theta = x(1);
        n = x(2);

        f = mu * theta^(1 - alpha);

        % Wage
        w = eta * (Z + c * theta) + (1 - eta) * b;

        % Job Value
        J = (Z - w) / A;

        % Free entry condition
        eq1 = c - beta * mu * theta^(-alpha) * J;

        % Employment steady-state condition
        eq2 = n - f / (delta + f);

        F = [eq1; eq2];
    end
end