diff --git a/andes/system.py b/andes/system.py
index 973b17fbd..379d84b02 100644
--- a/andes/system.py
+++ b/andes/system.py
@@ -13,6 +13,7 @@
 
 import configparser
 import importlib
+import importlib.util
 import inspect
 import logging
 import os
@@ -1248,7 +1249,6 @@ def connectivity(self, info=True):
         diag = list(matrix(spmatrix(u, to, os, (n, 1), 'd') +
                            spmatrix(u, fr, os, (n, 1), 'd')))
 
-        nib = self.Bus.n_islanded_buses = diag.count(0)
         for idx in range(n):
             if diag[idx] == 0:
                 self.Bus.islanded_buses.append(idx)
@@ -1265,12 +1265,20 @@ def connectivity(self, info=True):
                         'd')
         temp = sparse(temp)  # need to drop allocated zero values
 
-        cons = temp[0, :]
-        nelm = len(cons.J)
+        # Translated find_islanded_areas from goderya.jl into Python
         conn = spmatrix([], [], [], (1, n), 'd')
-        enum = idx = islands = 0
+        island_sets = []
+        starting_bus = 0
+        visit_idx = 0
 
         while True:
+            if starting_bus in self.Bus.islanded_buses:
+                starting_bus += 1
+                continue
+
+            cons = temp[starting_bus, :]
+            nelm = len(cons.J)
+
             while True:
                 cons = cons * temp
                 cons = sparse(cons)  # remove zero values
@@ -1279,31 +1287,24 @@ def connectivity(self, info=True):
                     break
                 nelm = new_nelm
 
-            # started with an islanded bus
-            if len(conn.J) == 0:
-                enum += 1
-            # all buses are interconnected
-            elif len(cons.J) == n:
-                break
-
-            self.Bus.island_sets.append(list(cons.J))
+            island_sets.append(list(cons.J))
             conn += cons
-            islands += 1
-            nconn = len(conn.J)
-            if nconn >= (n - nib):
-                self.Bus.island_sets = [i for i in self.Bus.island_sets if len(i) > 0]
+
+            if len(conn.J) >= (n - len(self.Bus.islanded_buses)):
                 break
 
-            for element in conn.J[idx:]:
-                if not diag[idx]:
-                    enum += 1  # skip islanded buses
-                if element <= enum:
-                    idx += 1
-                    enum += 1
+            # Increment `starting_bus` until it's not in `conn.J` and
+            # `self.Bus.islanded_buses`
+            for i in range(visit_idx, self.Bus.n):
+                if i in conn.J or i in self.Bus.islanded_buses:
+                    i += 1
                 else:
+                    visit_idx = i
                     break
 
-            cons = temp[enum, :]
+            starting_bus = visit_idx
+
+        self.Bus.island_sets = island_sets
 
         # --- check if all areas have a slack generator ---
         if len(self.Bus.island_sets) > 0: