diff --git a/deblur/deblur_2d.py b/deblur/deblur_2d.py
index 0676478..fd5da19 100644
--- a/deblur/deblur_2d.py
+++ b/deblur/deblur_2d.py
@@ -1,5 +1,7 @@
 import numpy as np
 from scipy.signal import convolve2d
+from scipy.sparse import lil_matrix
+from scipy.sparse.linalg import spsolve
 import cv2
 import matplotlib.pyplot as plt
 
@@ -25,6 +27,11 @@ print(kernel)
 # Perform 2D convolution (blurring)
 blurred = convolve2d(image, kernel, mode="same", boundary="fill", fillvalue=0)
 
+h, w = image.shape
+kh, kw = kernel.shape
+pad_h, pad_w = kh // 2, kw // 2
+
+
 show(image)
 show(blurred)
 
@@ -33,35 +40,28 @@ print("\nBlurred image:\n", blurred)
 
 print("\nBuilding linear system for deconvolution...")
 
-# Image size
-h, w = image.shape
-kh, kw = kernel.shape
-pad_h, pad_w = kh // 2, kw // 2
+# Step 2: Build sparse matrix A
+N = h * w
+A = lil_matrix((N, N), dtype=np.float32)
+b = blurred.flatten()
 
-# Build matrix A and vector b for Ax = b
-A = []
-b = []
+def index(y, x):
+    return y * w + x
 
 for y in range(h):
     for x in range(w):
-        row = np.zeros((h, w), dtype=np.float32)
-
+        row_idx = index(y, x)
         for ky in range(kh):
             for kx in range(kw):
                 iy = y + ky - pad_h
                 ix = x + kx - pad_w
                 if 0 <= iy < h and 0 <= ix < w:
-                    row[iy, ix] += kernel[ky, kx]
+                    col_idx = index(iy, ix)
+                    A[row_idx, col_idx] += kernel[ky, kx]
 
-        A.append(row.flatten())
-        b.append(blurred[y, x])
-
-A = np.array(A)
-b = np.array(b)
-
-# Solve for the deblurred image
-deblurred_flat = np.linalg.solve(A, b)
-deblurred = deblurred_flat.reshape((h, w))
+# Step 3: Solve the sparse system A * x = b
+x = spsolve(A.tocsr(), b)
+deblurred = x.reshape((h, w))
 
 print("\nDeblurred image:\n", np.round(deblurred, 2))