diff --git a/examples/demo.py b/examples/demo.py
new file mode 100644
index 0000000..6e38592
--- /dev/null
+++ b/examples/demo.py
@@ -0,0 +1,60 @@
+from numethods import *
+
+A = Matrix([[2, -1], [1, 2], [1, 1]])
+b = Vector([1, 2, 3])
+
+# Factorization
+qr = QRHouseholder(A)
+Q, R = qr.Q, qr.R
+print("Q =", Q)
+print("R =", R)
+
+qrm = QRModifiedGramSchmidt(A)
+Qm, Rm = qrm.Q, qrm.R
+print("Qm =", Qm)
+print("Rm =", Rm)
+print("Q^T Q =", Q.T @ Q)
+print("Qm^T Qm =", Qm.T @ Qm)
+print("A=Qm Rm =", Qm @ Rm)
+
+# Solve Ax = b (least squares, since A is tall)
+x_ls = LeastSquaresSolver(A, b).solve()
+print("Least squares solution:", x_ls)
+
+
+def demo_linear_solvers():
+    A = Matrix([[4, -1, 0], [-1, 4, -1], [0, -1, 3]])
+    b = Vector([15, 10, 10])
+
+    print("LU:", LUDecomposition(A).solve(b))
+    print("Gauss-Jordan:", GaussJordan(A).solve(b))
+    print("Cholesky:", Cholesky(A).solve(b))
+    print("Jacobi:", Jacobi(A, b, tol=1e-12).solve())
+    print("Gauss-Seidel:", GaussSeidel(A, b, tol=1e-12).solve())
+
+
+def demo_roots():
+    f = lambda x: x**3 - x - 2
+    df = lambda x: 3 * x**2 - 1
+    print("Bisection:", Bisection(f, 1, 2).solve())
+    print("Secant:", Secant(f, 1.0, 2.0).solve())
+    print("Newton root:", NewtonRoot(f, df, 1.5).solve())
+    # A simple contraction for demonstration; not general-purpose
+    g = lambda x: (x + 2 / x**2) / 2
+    print("Fixed point (demo):", FixedPoint(g, 1.5).solve())
+
+
+def demo_interpolation():
+    x = [0, 1, 2, 3]
+    y = [1, 2, 0, 5]
+    newt = NewtonInterpolation(x, y)
+    lagr = LagrangeInterpolation(x, y)
+    t = 1.5
+    print("Newton interpolation at", t, "=", newt.evaluate(t))
+    print("Lagrange interpolation at", t, "=", lagr.evaluate(t))
+
+
+if __name__ == "__main__":
+    demo_linear_solvers()
+    demo_roots()
+    demo_interpolation()