function [sigma, tangentModulus] = bilinMonotonic(E, fy, ks, strain)
    % Trial stress calculation
    sigma = E * strain;
    tangentModulus = E;

    % Plastic correction if stress exceeds yield stress
    if abs(sigma) > fy
        sigma = sign(sigma) * (fy + (abs(strain) - fy / E) * ks * E);
        tangentModulus = ks * E;
    end
end