diff --git a/Project/ApplyRobinBC_mult_rot.m b/Project/ApplyRobinBC_mult_rot.m
index 3a9d8ba..39b5587 100644
--- a/Project/ApplyRobinBC_mult_rot.m
+++ b/Project/ApplyRobinBC_mult_rot.m
@@ -1,36 +1,37 @@
 function [K, F] = ApplyRobinBC_mult_rot(model, K, F, alpha, u_out)
-mesh = model.Mesh;
+mesh  = model.Mesh;
 nodes = mesh.Nodes;
-elements = mesh.Elements;
 
-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,:);
+robinEdges = [4, 5, 6, 11];
 
-for k = 1:size(boundaryEdges,1)
-    i = boundaryEdges(k,1);
-    j = boundaryEdges(k,2);
+for eID = robinEdges
 
-    ri = nodes(1,i);
-    rj = nodes(1,j);
-    
-    %find the edges on the axis -> homogeneous Neumann BC, no Robin contribution
-    if ri == 0 && rj == 0
-        continue;
+    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
-
-    xi = nodes(:,i);
-    xj = nodes(:,j);
-
-    L = norm(xi - xj);
-    rbar = 0.5 * (xi(1) + xj(1));   % r at edge midpoint
-
-    Ke = rbar * alpha * L / 6 * [2 1; 1 2];
-    Fe = rbar * alpha * u_out * L / 2 * [1; 1];
-
-    K([i j],[i j]) = K([i j],[i j]) + Ke;
-    F([i j]) = F([i j]) + Fe;
 end
 end