diff --git a/examples/poly_eval/README.md b/examples/poly_eval/README.md
new file mode 100644
index 000000000..f4162db0b
--- /dev/null
+++ b/examples/poly_eval/README.md
@@ -0,0 +1,17 @@
+# Example - poly_eval 📘
+
+This example demonstrates the usage of the V Scientific Library for evaluating polynomial at given
+value of x.
+
+## Instructions
+
+1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io).
+2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)!
+3. Navigate to this directory.
+4. Run the example using the following command:
+
+```sh
+v run main.v
+```
+
+Enjoy exploring the capabilities of the V Scientific Library! 😊
diff --git a/examples/poly_eval/main.v b/examples/poly_eval/main.v
new file mode 100644
index 000000000..a82a51205
--- /dev/null
+++ b/examples/poly_eval/main.v
@@ -0,0 +1,16 @@
+module main
+
+import vsl.poly
+
+fn main() {
+	// represent the polynomial 2 * x^2 + 5 * x + 4
+	// with x = 4, result is 2 * 4^2 + 5 * 4 + 4 = 56
+	coef := [4.0, 5, 2]
+	result := poly.eval(coef, 4)
+	println('Evaluated value: ${result}')
+
+	// base polynomial: 2 * x^2 + 5 * x + 4 = 56 with x = 4
+	// 1st derivative: 4 * x + 5 = 21 with x = 4
+	result_derivs := poly.eval_derivs(coef, 4, 2)
+	println('Evaluated values: ${result_derivs}')
+}
diff --git a/examples/poly_matrix/README.md b/examples/poly_matrix/README.md
new file mode 100644
index 000000000..76b0e48d9
--- /dev/null
+++ b/examples/poly_matrix/README.md
@@ -0,0 +1,16 @@
+# Example - poly_matrix 📘
+
+This example demonstrates the usage of the V Scientific Library for constructing companion matrix.
+
+## Instructions
+
+1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io).
+2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)!
+3. Navigate to this directory.
+4. Run the example using the following command:
+
+```sh
+v run main.v
+```
+
+Enjoy exploring the capabilities of the V Scientific Library! 😊
diff --git a/examples/poly_matrix/main.v b/examples/poly_matrix/main.v
new file mode 100644
index 000000000..23a966e14
--- /dev/null
+++ b/examples/poly_matrix/main.v
@@ -0,0 +1,27 @@
+module main
+
+import vsl.poly
+
+fn main() {
+	// Polynomial coefficients for p(x) = 2x^3 -4x^2 + 3x - 5
+	coef := [2.0, -4.0, 3.0, -5.0]
+
+	// For the polynomial p(x) = a_n x^n + a_{n-1} x^{n-1} + ... + a_1 x + a_0
+	// The companion matrix C will be:
+	// [ 0    0    0   -a_0/a_n ]
+	// [ 1    0    0   -a_1/a_n ]
+	// [ 0    1    0   -a_2/a_n ]
+	// [ 0    0    1   -a_3/a_n ]
+	comp_matrix := poly.companion_matrix(coef)
+	println('Companion matrix:')
+	for row in comp_matrix {
+		println(row)
+	}
+
+	// Balancing matrix if needed
+	balanced_matrix := poly.balance_companion_matrix(comp_matrix)
+	println('Balanced companion matrix:')
+	for row in balanced_matrix {
+		println(row)
+	}
+}
diff --git a/examples/poly_operations/README.md b/examples/poly_operations/README.md
new file mode 100644
index 000000000..79e818955
--- /dev/null
+++ b/examples/poly_operations/README.md
@@ -0,0 +1,17 @@
+# Example - poly_eval 📘
+
+This example demonstrates the usage of the V Scientific Library for additioning, substracting and
+multiplying polynomials.
+
+## Instructions
+
+1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io).
+2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)!
+3. Navigate to this directory.
+4. Run the example using the following command:
+
+```sh
+v run main.v
+```
+
+Enjoy exploring the capabilities of the V Scientific Library! 😊
diff --git a/examples/poly_operations/main.v b/examples/poly_operations/main.v
new file mode 100644
index 000000000..8985948d5
--- /dev/null
+++ b/examples/poly_operations/main.v
@@ -0,0 +1,40 @@
+module main
+
+import vsl.poly
+
+fn main() {
+	// Addition
+	// Degree is not modified unless highest coefficients cancel each other out
+	poly_1 := [1.0, 12, 3] // 1 + 12x + 3x^2
+	poly_2 := [3.0, 2, 7] // 3 + 2x + 7x^2
+	result_add := poly.add(poly_1, poly_2)
+	println('Addition result: ${result_add}') // Expected: [4.0, 14.0, 10.0] (4 + 14x + 10x^2)
+
+	// Subtraction
+	// Degree is not modified unless highest coefficients cancel each other out
+	result_sub := poly.subtract(poly_1, poly_2)
+	println('Subtraction result: ${result_sub}') // Expected: [-2.0, 10.0, -4.0] (-2 + 10x - 4x^2)
+
+	// Multiplication
+	// with given degree n and m for poly_1 and poly_2
+	// resulting polynomial will be of degree n + m
+	result_mult := poly.multiply(poly_1, poly_2)
+	println('Multplication result: ${result_mult}') // Expected: [3.0, 38.0, 40.0, 90.0, 21.0] (3 + 38x + 400x^2 + 90x^3 + 21x^4)
+
+	// Division
+	// Result includes the quotient and the remainder
+	// To get the real remainder, divide it by the divisor.
+
+	// OLD WAY
+	// poly_dividend := [2.0, -4.0, -4.0, 1.0] // 2 - 4x - 4x^2 + x^3
+	// poly_divisor := [-2.0, 1.0] // -2 + x
+
+	// REVERSED WAY
+	poly_dividend := [2.0, -4.0, -4.0, 1.0] // 2 - 4x - 4x^2 + x^3
+	poly_divisor := [-2.0, 1.0] // -2 + x
+
+	quotient, remainder := poly.divide(poly_dividend, poly_divisor)
+	println('Division quotient: ${quotient}') // Expected quotient: [-8.0, -2.0, 1.0] (-8 - 2x + x^2)
+	println('Division remainder: ${remainder}') // Expected remainder: [-14.0]
+	// Real remainder: -14 / (-2 + x)
+}
diff --git a/poly/poly.v b/poly/poly.v
index acf36c397..fdaaa70f0 100644
--- a/poly/poly.v
+++ b/poly/poly.v
@@ -15,9 +15,11 @@ pub fn eval(c []f64, x f64) f64 {
 		errors.vsl_panic('coeficients can not be empty', .efailed)
 	}
 	len := c.len
-	mut ans := c[len - 1]
-	for e in c[..len - 1] {
-		ans = e + x * ans
+	mut ans := 0.0
+	mut i := len - 1
+	for i >= 0 {
+		ans = c[i] + x * ans
+		i--
 	}
 	return ans
 }
@@ -144,13 +146,13 @@ fn sorted_3_(x_ f64, y_ f64, z_ f64) (f64, f64, f64) {
 	mut y := y_
 	mut z := z_
 	if x > y {
-		y, x = swap_(x, y)
+		x, y = swap_(x, y)
 	}
-	if y > z {
-		z, y = swap_(y, z)
+	if x > z {
+		x, z = swap_(x, z)
 	}
-	if x > y {
-		y, x = swap_(x, y)
+	if y > z {
+		y, z = swap_(y, z)
 	}
 	return x, y, z
 }
@@ -213,17 +215,17 @@ pub fn balance_companion_matrix(cm [][]f64) [][]f64 {
 			if col_norm == 0.0 || row_norm == 0.0 {
 				continue
 			}
-			mut g := row_norm / poly.radix
+			mut g := row_norm / radix
 			mut f := 1.0
 			s := col_norm + row_norm
 			for col_norm < g {
-				f *= poly.radix
-				col_norm *= poly.radix2
+				f *= radix
+				col_norm *= radix2
 			}
-			g = row_norm * poly.radix
+			g = row_norm * radix
 			for col_norm > g {
-				f /= poly.radix
-				col_norm /= poly.radix2
+				f /= radix
+				col_norm /= radix2
 			}
 			if (row_norm + col_norm) < 0.95 * s * f {
 				not_converged = true
@@ -288,25 +290,34 @@ pub fn multiply(a []f64, b []f64) []f64 {
 // Output: (q, r) where q is the quotient and r is the remainder
 // such that a(x) = b(x) * q(x) + r(x) and degree(r) < degree(b)
 pub fn divide(a []f64, b []f64) ([]f64, []f64) {
-	mut quotient := []f64{}
+	if b.len == 0 {
+		panic('divisor cannot be an empty polynomial')
+	}
+	if a.len == 0 {
+		return []f64{len: 0}, []f64{len: 0}
+	}
+
+	mut quotient := []f64{len: a.len - b.len + 1, init: 0.0}
 	mut remainder := a.clone()
-	b_lead_coef := b[0]
+
+	b_degree := b.len - 1
+	b_lead_coeff := b[b_degree]
 
 	for remainder.len >= b.len {
-		lead_coef := remainder[0] / b_lead_coef
-		quotient << lead_coef
+		remainder_degree := remainder.len - 1
+		lead_coeff := remainder[remainder_degree]
+
+		quotient_term := lead_coeff / b_lead_coeff
+		quotient_idx := remainder_degree - b_degree
+		quotient[quotient_idx] = quotient_term
+
 		for i in 0 .. b.len {
-			remainder[i] -= lead_coef * b[i]
-		}
-		remainder = unsafe { remainder[1..] }
-		for remainder.len > 0 && math.abs(remainder[0]) < 1e-10 {
-			remainder = unsafe { remainder[1..] }
+			remainder[quotient_idx + i] -= quotient_term * b[i]
 		}
-	}
 
-	if remainder.len == 0 {
-		remainder = []f64{}
+		for remainder.len > 0 && remainder[remainder.len - 1] == 0.0 {
+			remainder = remainder[0..remainder.len - 1].clone()
+		}
 	}
-
 	return quotient, remainder
 }
diff --git a/poly/poly_test.v b/poly/poly_test.v
index b9c25cdb2..82571b6eb 100644
--- a/poly/poly_test.v
+++ b/poly/poly_test.v
@@ -1,14 +1,24 @@
 module poly
 
-import math
-
 fn test_eval() {
-	// ans = 2
-	// ans = 4.0 + 4 * 2 = 12
-	// ans = 5 + 4 * 12 = 53
 	x := 4
-	cof := [4.0, 5, 2]
-	assert eval(cof, 4) == 53
+	coef := [4.0, 5, 2]
+	assert eval(coef, 4) == 56
+}
+
+fn test_eval_derivs() {
+	coef := [4.0, 3.0, 2.0]
+	x := 1.0
+	lenres := 3
+	assert eval_derivs(coef, x, lenres) == [9.0, 7.0, 4.0]
+}
+
+fn test_solve_quadratic() {
+	a := 1.0
+	b := -3.0
+	c := 2.0
+	roots := solve_quadratic(a, b, c)
+	assert roots == [1.0, 2.0]
 }
 
 fn test_swap() {
@@ -18,11 +28,37 @@ fn test_swap() {
 	assert a == 202.0 && b == 101.0
 }
 
+fn test_sorted_3_() {
+	a := 5.0
+	b := 7.0
+	c := -8.0
+	x, y, z := sorted_3_(a, b, c)
+	assert x == -8.0 && y == 5.0 && z == 7.0
+}
+
+fn test_companion_matrix() {
+	coef := [2.0, -4.0, 3.0, -5.0]
+	assert companion_matrix(coef) == [
+		[0.0, 0.0, 0.4],
+		[1.0, 0.0, -0.8],
+		[0.0, 1.0, 0.6],
+	]
+}
+
+fn test_balance_companion_matrix() {
+	coef := [2.0, -4.0, 3.0, -5.0]
+	matrix := companion_matrix(coef)
+	assert balance_companion_matrix(matrix) == [
+		[0.0, 0.0, 0.8],
+		[0.5, 0.0, -0.8],
+		[0.0, 1.0, 0.6],
+	]
+}
+
 fn test_add() {
 	a := [1.0, 2.0, 3.0]
 	b := [6.0, 20.0, -10.0]
 	result := add(a, b)
-	println('Add result: ${result}')
 	assert result == [7.0, 22.0, -7.0]
 }
 
@@ -30,7 +66,6 @@ fn test_subtract() {
 	a := [6.0, 1532.0, -4.0]
 	b := [1.0, -1.0, -5.0]
 	result := subtract(a, b)
-	println('Subtract result: ${result}')
 	assert result == [5.0, 1533.0, 1.0]
 }
 
@@ -39,17 +74,19 @@ fn test_multiply() {
 	a := [2.0, 3.0, 4.0]
 	b := [0.0, -3.0, 2.0]
 	result := multiply(a, b)
-	println('Multiply result: ${result}')
 	assert result == [0.0, -6.0, -5.0, -6.0, 8.0]
 }
 
 fn test_divide() {
-	// (x^2 + 2x + 1) / (x + 1) = (x + 1)
-	a := [1.0, 2.0, 1.0]
-	b := [1.0, 1.0]
-	quotient, remainder := divide(a, b)
-	println('Divide quotient: ${quotient}')
-	println('Divide remainder: ${remainder}')
-	assert quotient == [1.0, 1.0]
-	assert remainder == [] // The empty set indicates that two polynomials divide each other exactly (without remainder).
+	a := [0.0, 0.0, 1.0, 1.0]
+	b := [-2.0, 1.0, 1]
+	mut quotient, mut remainder := divide(a, b)
+	assert quotient == [0.0, 1.0]
+	assert remainder == [0.0, 2.0]
+
+	c := [5.0, -11.0, -7.0, 4.0]
+	d := [5.0, 4.0]
+	quotient, remainder = divide(c, d)
+	assert quotient == [1.0, -3.0, 1]
+	assert remainder == []
 }