function [K, F] = ApplyRobinBC_mult(model, K, F, alpha, u_out)
% ApplyRobinBC_mult
% Implements Robin BC:
%   lambda * du/dn + alpha*u = alpha*u_out
%
% alpha   : heat transfer coefficient
% u_out   : outside temperature

mesh = model.Mesh;
nodes = mesh.Nodes;
elements = mesh.Elements;

% --- Find boundary edges ---
edgesAll = [elements([1 2],:), elements([2 3],:), elements([3 1],:)];
edgesSorted = sort(edgesAll,1);
[edgesUnique,~,ic] = unique(edgesSorted','rows');
counts = accumarray(ic,1);
boundaryEdges = edgesUnique(counts==1,:);   % Nx2 array

% --- Loop over Robin boundary edges ---
for k = 1:size(boundaryEdges,1)
    i = boundaryEdges(k,1);
    j = boundaryEdges(k,2);

    ri = nodes(1,i);
    rj = nodes(1,j);
    
    if ri == 0 && rj == 0
    % Edge lies on symmetry axis r = 0
    % -> homogeneous Neumann BC, no Robin contribution
        continue;
    end

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

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

    % Robin boundary element matrices
    Ke = alpha * L / 6 * [2 1; 1 2];
    Fe = 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