function [K, F] = ApplyRobinBC_mult_rot(model, K, F, alpha, u_out)
mesh  = model.Mesh;
nodes = mesh.Nodes;

robinEdges = [4, 5, 6, 11];

for eID = robinEdges

    edgeNodes = findNodes(mesh,'region','Edge',eID);
    for k = 1:length(edgeNodes)-1
        i = edgeNodes(k);
        j = edgeNodes(k+1);

        xi = nodes(:,i);
        xj = nodes(:,j);

        % Edge length
        L = norm(xi - xj);

        % Axisymmetric radius at midpoint
        rbar = 0.5 * (xi(1) + xj(1));

        % Skip axis (safety, should not occur here)
        if rbar == 0
            continue;
        end

        % Robin element matrices (linear edge)
        Ke = rbar * alpha * L / 6 * [2 1; 1 2];
        Fe = rbar * alpha * u_out * L / 2 * [1; 1];

        % Assemble
        K([i j],[i j]) = K([i j],[i j]) + Ke;
        F([i j])      = F([i j])      + Fe;
    end
end
end