diff --git a/lab3_page7.py b/lab3_page7.py
new file mode 100644
index 0000000..3a1f6c8
--- /dev/null
+++ b/lab3_page7.py
@@ -0,0 +1,139 @@
+# ==== Lab #3: Convolution & Sinusoids ====
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.integrate import solve_ivp
+
+# Components (same as Lab #2)
+R1, C1 = 10e3, 1e-6
+R2, C2 = 4.7e3, 1e-6
+tau1, tau2 = R1*C1, R2*C2
+R4, R5 = 1e3, 4.7e3
+G = 1 + R5/R4
+
+# Time grid
+tau_min = min(tau1, tau2)
+tau_max = max(tau1, tau2)
+dt = 0.0001
+t = np.arange(0, 5 + dt/2, dt)
+N = t.size
+
+# Impulse responses (causal)
+h1 = (1/tau1) * np.exp(-t/tau1) # u(t) implicit (0 for t<0)
+h2 = (1/tau2) * np.exp(-t/tau2)
+
+# Overall impulse response for cascaded stages (discrete conv scaled by dt)
+h = G * dt * np.convolve(h1, h2, mode='full')
+
+# Helper: numeric convolution with step dt (Riemann sums)
+def conv_num(a, b, dt):
+    """Full-length discrete convolution approximating integrals via Riemann sums with step dt."""
+    return np.convolve(a, b, mode='full') * dt
+
+# ======================
+# 1) Single sinusoids
+# ======================
+A = 1.0
+freqs = [1, 5, 10, 50]
+xmax = 2.0
+
+plt.figure(1)
+plt.clf()
+for k, f in enumerate(freqs, start=1):
+    w = 2*np.pi*f
+    x = A * np.sin(w*t) * (t >= 0)
+    y_full = conv_num(x, h, dt)
+    y = y_full[:N] # keep same length as x
+
+# ODE: tau1*tau2*y'' + (tau1+tau2)*y' + y = G*x(t)
+def rhs(ti, z):
+    x_t = A*np.sin(w*ti) * (ti >= 0.0)
+    return [
+        z[1],
+        (G*x_t - (tau1 + tau2)*z[1] - z[0]) / (tau1*tau2)
+    ]
+z0 = [0.0, 0.0] # zero initial conditions
+sol = solve_ivp(rhs, (t[0], t[-1]), z0, t_eval=t, rtol=1e-9, atol=1e-12)
+y_ode = sol.y[0]
+
+ax = plt.subplot(len(freqs), 1, k)
+ax.plot(t, x, label='input')
+ax.plot(t, y, label='output (h*x)')
+ax.plot(t, y_ode, '--', label='output (ODE)', color='r')
+ax.grid(True)
+ax.set_xlabel('t (s)')
+ax.set_ylabel('V')
+ax.set_xlim([0, xmax])
+ax.set_title(f'Single sinusoid: f = {f:g} Hz')
+ax.legend()
+xmax = xmax/2
+plt.tight_layout()
+
+# =========================================
+# 2) Sum of sinusoids
+# =========================================
+A1, f1 = 1.0, 5.0
+A2, f2 = 0.5, 20.0
+x1 = A1 * np.sin(2*np.pi*f1*t) * (t >= 0)
+x2 = A2 * np.sin(2*np.pi*f2*t) * (t >= 0)
+alpha, beta = 1.5, 2.0
+x = alpha*x1 + beta*x2
+Nx = x.size
+
+y = conv_num(x, h, dt)[:Nx]
+y1 = conv_num(x1, h, dt)[:Nx]
+y2 = conv_num(x2, h, dt)[:Nx]
+
+plt.figure(2)
+plt.clf()
+plt.plot(t, y, label='y = h*(x1+x2)')
+plt.plot(t, alpha*y1 + beta*y2, '--', label='alpha*y1 + beta*y2')
+plt.grid(True)
+plt.xlabel('t (s)')
+plt.ylabel('V')
+plt.xlim([0, 1])
+plt.legend()
+plt.title('Linearity check: y[u1+u2] vs y[u1]+y[u2]')
+
+# ==================================
+# 3) Time invariance (delay \Delta)
+# ==================================
+Delta = 0.1
+nD = int(round(Delta/dt))
+x_del = np.concatenate([np.zeros(nD), x[:-nD]]) if nD < Nx else np.zeros_like(x)
+y_del = conv_num(x_del, h, dt)[:x_del.size]
+
+plt.figure(3)
+plt.clf()
+plt.plot(t, y, label='y(t)')
+plt.plot(t, np.concatenate([np.zeros(nD), y[:-nD]]), label=r'$y(t-\Delta)$')
+plt.plot(t, y_del, '--', label=r'$h*u(t-\Delta)$')
+plt.grid(True)
+plt.xlabel('t (s)')
+plt.ylabel('V')
+plt.xlim([0, 1])
+plt.legend()
+plt.title('Time invariance check')
+plt.tight_layout()
+
+
+# =======================================
+# 4) Commutativity of stages (h1*h2)
+# =======================================
+h12 = G * dt * np.convolve(h1, h2, mode='full')
+h21 = G * dt * np.convolve(h2, h1, mode='full')
+rmse = np.sqrt(np.mean((h12 - h21)**2))
+
+plt.figure(4)
+plt.clf()
+plt.plot(t, h12[:N], label=r'$h_1 \ast h_2$')
+plt.plot(t, h21[:N], '--', label=r'$h_2 \ast h_1$')
+plt.grid(True)
+plt.xlim([0, 0.1])
+plt.xlabel('t (s)')
+plt.ylabel('V')
+plt.legend()
+plt.title(rf'Commutativity - RMSE of $ e = h_{{12}} - h_{{21}}: \sqrt{{\mathrm{{mean}}(e^2)}} = {rmse:.3e} $')
+
+plt.tight_layout()
+plt.show()