diff --git a/project_euler/problem_061/__init__.py b/project_euler/problem_061/__init__.py
new file mode 100644
index 000000000000..e69de29bb2d1
diff --git a/project_euler/problem_061/sol1.py b/project_euler/problem_061/sol1.py
new file mode 100644
index 000000000000..c4146f4f0364
--- /dev/null
+++ b/project_euler/problem_061/sol1.py
@@ -0,0 +1,133 @@
+"""
+Project Euler Problem 61: https://projecteuler.net/problem=61
+
+Triangle, square, pentagonal, hexagonal, heptagonal, and octagonal numbers are a
+ll figurate (polygonal) numbers and are generated by the following formulae:
+
+Triangle    P(3,n)=n(n+1)/2         1, 3, 6, 10, 15, ...
+Square 		P(4,n)=n^2 		        1, 4, 9, 16, 25, ...
+Pentagonal 	P(5,n)=n(3n-1)/2 		1, 5, 12, 22, 35, ...
+Hexagonal 	P(6,n)=n(2n-1) 		    1, 6, 15, 28, 45, ...
+Heptagonal 	P(7,n)=n(5n-3)/2 		1, 7, 18, 34, 55, ...
+Octagonal 	P(8,n)=n(3n-2) 		    1, 8, 21, 40, 65, ...
+
+The ordered set of three 4-digit numbers: 8128, 2882, 8281, has
+three interesting properties.
+
+    1. The set is cyclic, in that the last two digits of each number is the first two
+        digits of the next number (including the last number with the first).
+    2. Each polygonal type: triangle (P(3,127) = 8128), square (P(4,91) = 8281), and
+        pentagonal (P(5,44) = 2882),is represented by a different number in the set.
+    3. This is the only set of 4-digit numbers with this property.
+
+Find the sum of the only ordered set of six cyclic 4-digit numbers for which each
+polygonal type: triangle, square, pentagonal, hexagonal, heptagonal, and octagonal,
+is represented by a different number in the set.
+"""
+
+from itertools import permutations
+
+
+def solution() -> int:
+    """
+    For this task is good to know some basics of combinatorial analysis, and
+    know how to solve a quadratic equation.
+
+    On the first sight we can see that number of possibilities is quite big. But
+    from the description of the task, we can get to know that:
+        1. We are only looking for 4 digits numbers
+        2. All numbers are naturals numbers.
+
+    By knowing that we can try to limit range of numbers that
+    we will be looking for specific formulae.
+
+    We need to solve some quadratic equation to get that range:
+        1. Triangle:
+            Lower: n(n+1)/2 = 1000 // 45
+
+            ######
+            We get two solutions one positive, one negative. From previous observation
+            we take positive one, and round it down or up, so
+            we get maximum 4 digit result. We doing that for rest of formulae.
+            ######
+
+            Upper: n(n+1)/2 = 10000 // 140
+
+        2. Square:
+            Lower: n^2 = 1000 // 32
+            Upper: n^2 = 9999 // 99
+        3. Pentagonal:
+            Lower: n(3n-1)/2 = 1000 // 26
+            Upper: n(3n-1)/2 = 9999 // 81
+        4. Hexagonal:
+            Lower: n(2n-1) = 1000 // 23
+            Upper: n(2n-1) = 9999 // 70
+        5. Heptagonal:
+            Lower: n(5n-3)/2 = 1000 // 21
+            Upper: n(5n-3)/2 = 9999 // 63
+        6. Octagonal
+            Lower: n(3n-2) = 1000 // 19
+            Upper: n(3n-2) = 9999 // 59
+
+    Then we just need to check if the last two digits of each number
+    is the first two digits of the next number, and if they are not the same.
+
+    I created a function is_cyclic to check if it the last two digits of number is
+    the first two digits of the next number.
+
+    Then program iterates through all permutations of polygonal types. (itertools).
+    For each permutation, it iterates through the corresponding lists of
+    polygonal numbers stored in the polygonals dictionary.
+
+    At the end if 6 polygonal numbers form a cyclic set, it returns the sum of them.
+    """
+
+    triangle = [int(x * (x + 1) * 0.5) for x in range(45, 141)]
+
+    square = [int(x * x) for x in range(32, 100)]
+
+    pentagonal = [int(x * (3 * x - 1) * 0.5) for x in range(26, 82)]
+
+    hexagonal = [int(x * (2 * x - 1)) for x in range(23, 71)]
+
+    heptagonal = [int(x * (5 * x - 3) * 0.5) for x in range(21, 64)]
+
+    octagonal = [int(x * (3 * x - 2)) for x in range(19, 60)]
+
+    polygonals = {
+        3: triangle,
+        4: square,
+        5: pentagonal,
+        6: hexagonal,
+        7: heptagonal,
+        8: octagonal,
+    }
+
+    for perm in permutations(range(3, 9)):
+        for t in polygonals[perm[0]]:
+            for s in polygonals[perm[1]]:
+                if not is_cyclic(t, s):
+                    continue
+                for p in polygonals[perm[2]]:
+                    if not is_cyclic(s, p):
+                        continue
+                    for hx in polygonals[perm[3]]:
+                        if not is_cyclic(p, hx):
+                            continue
+                        for hp in polygonals[perm[4]]:
+                            if not is_cyclic(hx, hp):
+                                continue
+                            for o in polygonals[perm[5]]:
+                                if not is_cyclic(hp, o) or not is_cyclic(o, t):
+                                    continue
+                                return sum([t, s, p, hx, hp, o])
+
+    return 0
+
+
+def is_cyclic(a: int, b: int):
+    return str(a)[2:] == str(b)[:2]
+
+
+if __name__ == "__main__":
+    print(f"{solution() = }")