diff --git a/code/main.py b/code/main.py
index cb64fc1..50bbc72 100644
--- a/code/main.py
+++ b/code/main.py
@@ -12,6 +12,10 @@ def create_subplot():
    return plt.subplots(layout='constrained')[1]
 
 
+def plt_append(sp, x: list[float], y: list[float], label: str, format: str):
+   sp.plot(x, y, format, label=label)
+
+
 class NonLinear:
    bisect_exp = "x**2 * np.sin(x)"
    newton_exp = "np.sin(x) * np.sqrt(np.abs(x))"
@@ -38,10 +42,6 @@ class NonLinear:
           range, val_max) if val_max is not None else range.index(max(range))
       return range[index_l:index_r+1]
 
-   @staticmethod
-   def plt_append(sp, x: list[float], y: list[float], label: str, format: str):
-      sp.plot(x, y, format, label=label)
-
    @staticmethod
    def bisect(x, x_min, x_max):
       def f(x): return eval(NonLinear.bisect_exp)
@@ -62,11 +62,11 @@ class NonLinear:
       sol1 = NonLinear.bisect(x1, bounds[0], bounds[1])
       sol2 = NonLinear.bisect(x2, split_val, bounds[1])
 
-      NonLinear.plt_append(
+      plt_append(
           sp, x1, sol1[0], f"Исходные данные (y={NonLinear.bisect_exp})", "-b")
-      NonLinear.plt_append(
+      plt_append(
           sp, *(sol1[1]), f"bisect at [{bounds[0]},{bounds[1]}]", "or")
-      NonLinear.plt_append(
+      plt_append(
           sp, *(sol2[1]), f"bisect at [{split_val},{bounds[1]}]", "og")
 
       sp.set_title("scipy.optimize.bisect")
@@ -92,11 +92,11 @@ class NonLinear:
       sol1 = NonLinear.newton(x1, x0_1)
       sol2 = NonLinear.newton(x2, x0_2)
 
-      NonLinear.plt_append(
+      plt_append(
           sp, x1, sol1[0], f"Исходные данные (y={NonLinear.newton_exp})", "-b")
-      NonLinear.plt_append(
+      plt_append(
           sp, *(sol1[1]), f"newton at [{bounds[0]},{bounds[1]}]", "or")
-      NonLinear.plt_append(
+      plt_append(
           sp, *(sol2[1]), f"newton at [{split_l},{bounds[1]}]", "og")
 
       sp.set_title("scipy.optimize.newton")
@@ -113,8 +113,112 @@ class NonLinear:
       plt.show()
 
 
+class SLE:
+   gauss_data = ([[13, 2], [3, 4]], [1, 2])
+   invmatrix_data = ([[13, 2], [3, 4]], [1, 2])
+   tridiagonal_data = ([[4,  5,  6,  7, 8, 9],
+                       [2,  2,  2,  2, 2, 0]],
+                       [1, 2, 2, 3, 3, 3])
+
+   @staticmethod
+   def var_str(index):
+      return f"x{index+1}"
+
+   @staticmethod
+   def print_solution(data: list[float]):
+      print("  ", end='')
+      for i, val in enumerate(data[:-1]):
+         print(f"{SLE.var_str(i)} = {round(val,3)}, ", end='')
+      print(f"{SLE.var_str(len(data)-1)} = {round(data[-1],3)}")
+
+   @staticmethod
+   def print_data(data: tuple[list[list[float]], list[float]], tridiagonal: bool = False):
+      if tridiagonal:
+         new_data = []
+         new_len = len(data[0][0])
+         zipped = list(zip(*tuple(data[0])))
+         zipped[len(zipped)-1] = (zipped[len(zipped)-1][0],zipped[len(zipped)-2][1])
+         complement_to = new_len - len(zipped[0])
+         for i, val in enumerate(zipped):
+            zero_r = complement_to - i
+            if zero_r <= 0:
+               zero_r = 0
+            mid_val = list(reversed(val[1:])) + list(val)
+            mid_end = len(mid_val) if zero_r > 0 else len(
+                mid_val) + (complement_to - i)
+            mid_beg = len(mid_val) - (new_len - zero_r) if zero_r > 0 else 0
+            mid_beg = mid_beg if mid_beg >= 0 else 0
+            zero_l = new_len - (zero_r + (mid_end - mid_beg))
+            tmp = [0] * zero_l + \
+                mid_val[mid_beg:mid_end] + [0] * zero_r
+            new_data.append(tmp)
+         data = (new_data, data[1])
+      for i, val in enumerate(data[0]):
+         print("  ", end='')
+         for i_coef, coef in enumerate(val[:-1]):
+            if coef != 0:
+               print(f"({coef}{SLE.var_str(i_coef)}) + ", end='')
+            else:
+               print(f"  {coef}   + ",end='')
+         print(f"({val[-1]}{SLE.var_str(len(val)-1)})", end='')
+         print(f" = {data[1][i]}")
+
+   @staticmethod
+   def gauss(system: list[list[float]], b: list[float]):
+      lup = salg.lu_factor(system)
+      solution = salg.lu_solve(lup, b)
+      return solution
+
+   @staticmethod
+   def invmatrix(system: list[list[float]], b: list[float]):
+      m_inv = salg.inv(system)
+      solution = m_inv @ b
+      return solution
+
+   @staticmethod
+   def tridiagonal(system: list[list[float]], b: list[float]):
+      solution = salg.solveh_banded(system, b, lower=True)
+      return solution
+
+   @staticmethod
+   def print_gauss():
+      print("Gauss method (LU decomposition)")
+      print("  Input system:")
+      SLE.print_data(SLE.gauss_data)
+      print("  Solution:")
+      SLE.print_solution(SLE.gauss(*SLE.gauss_data))
+
+   @staticmethod
+   def print_invmatrix():
+      print("Inverted matrix method")
+      print("  Input system:")
+      SLE.print_data(SLE.invmatrix_data)
+      print("  Solution:")
+      SLE.print_solution(SLE.invmatrix(*SLE.invmatrix_data))
+
+   @staticmethod
+   def print_tridiagonal():
+      print("Tridiagonal matrix method (Thomas algorithm)")
+      print("  Input system:")
+      SLE.print_data(SLE.tridiagonal_data, True)
+      print("  Solution:")
+      SLE.print_solution(SLE.tridiagonal(*SLE.tridiagonal_data))
+
+   @staticmethod
+   def print(method="all"):
+      if method in ["gauss", "all"]:
+         SLE.print_gauss()
+      if method in ["invmatrix", "all"]:
+         SLE.print_invmatrix()
+      if method in ["banded", "all"]:
+         SLE.print_tridiagonal()
+
+class Approx:
+   function = "np.sin(x) * np.sqrt(np.abs(x))"
+
 def main():
-   NonLinear.plot()
+   # NonLinear.plot()
+   SLE.print()
 
 
 if __name__ == "__main__":