diff --git a/bin/abc_lab b/bin/abc_lab
index 2300648..1dcbf55 100755
Binary files a/bin/abc_lab and b/bin/abc_lab differ
diff --git a/include/utils/Grid1D.h b/include/utils/Grid1D.h
new file mode 100644
index 0000000..bd91615
--- /dev/null
+++ b/include/utils/Grid1D.h
@@ -0,0 +1,39 @@
+#ifndef _grid1d_n_
+#define _grid1d_n_
+
+#include "./utils/matrix.h"
+
+
+namespace utils{
+	//#######################################
+	//#            Grid1D TYPE              #
+	//#######################################
+	template <typename T>
+	struct Grid1D{
+		utils::Vector<T> grid;
+		utils::Vector<T> vertices;
+		utils::Vector<T> vertices_norm;
+
+	void create_vertices_norm(){
+		vertices_norm.fill(vertices.size()*2, 0);
+		uint64_t k = 0;
+		for (uint64_t i = 0; i < grid.size(); i++){
+			for (uint64_t j = 1; j <= 2; j++){
+				vertices_norm[k] = grid[i] - vertices[i+j];
+				k++;
+			}
+			
+			//vertices_norm[(i*2)+1] = grid[i] - vertices[(i*2)+1];
+		}
+		vertices_norm.print();
+	}
+};
+
+
+	typedef Grid1D<int> Gridi;
+	typedef Grid1D<float> Gridf;
+	typedef Grid1D<double> Gridd;
+}
+
+
+#endif // _grid1d_n_
\ No newline at end of file
diff --git a/include/utils/matrix.h b/include/utils/matrix.h
index 1fe32f0..7fb1abb 100644
--- a/include/utils/matrix.h
+++ b/include/utils/matrix.h
@@ -1,10 +1,6 @@
 #ifndef _matrix_n_
 #define _matrix_n_
 
-#include "iostream"
-#include <vector>
-//#include <random>
-
 #include "./utils/vector.h"
 
 
@@ -16,6 +12,8 @@ namespace utils{
 	struct Matrix{
 		utils::Vector<T> m;
 
+		T& operator[](uint64_t idx) { return m[idx]; } // Makes it able to do matr[1][1]
+		const T& operator[](uint64_t idx) const { return m[idx]; } // Makes it able to do matr[1][1]
 
 		using vector_type = typename decltype(std::declval<T>().v)::value_type;
 
@@ -26,8 +24,7 @@ namespace utils{
 
 	    Matrix() = default; // Default constructor
 
-		T& operator[](uint64_t idx) { return m[idx]; } // Makes it able to do matr[1][1]
-		const T& operator[](uint64_t idx) const { return m[idx]; } // Makes it able to do matr[1][1]
+
 
 
 		void fill(uint64_t rows, uint64_t cols, const vector_type num=0){
diff --git a/include/utils/utils.h b/include/utils/utils.h
index cfd11c5..bbef9ab 100644
--- a/include/utils/utils.h
+++ b/include/utils/utils.h
@@ -3,3 +3,4 @@
 
 #include "./utils/vector.h"
 #include "./utils/matrix.h"
+#include "./utils/Grid1D.h"
diff --git a/include/utils/vector.h b/include/utils/vector.h
index 80ad505..77a39ec 100644
--- a/include/utils/vector.h
+++ b/include/utils/vector.h
@@ -7,91 +7,196 @@
 
 
 namespace utils{
-	//#######################################
+	//######################################
 	//#            VECTOR TYPE             #
-	//#######################################
+	//# 	 Backed by std::vector<T>      #
+	//######################################
 	template<typename T>
-	struct Vector{
-
+	class Vector{
+public:
 		using value_type = T;
 		std::vector<T> v;
 
-	    Vector() = default;
+		// Default construtor utils::Vd vector;
+	    Vector() = default; 
 
-	    Vector(uint64_t size, T value = T()) {
+	    // Construtor utils::Vd vector(2, 3.15);
+	    // Prevent implicit conversions like Vector<double> v = 5;
+	    explicit Vector(uint64_t size, T value = T()) { 
 	        v.resize(size, value);
 	    }
 
 
-		void clear(){v.clear();}
-		void push_back(const T& val) { v.push_back(val); }
+	    //##########################################################
+		//#            VECTOR: --- basic properties ---            #
+		//##########################################################
+		// vector.clear();
+		void clear() noexcept {v.clear();}
+
+		// vector.push_back(2);
+		void push_back(const T& val) { v.push_back(val); } 
+
+		// vector[2] = 2;
 		T& operator[](uint64_t idx) { return v[idx]; }
+		// a = vector[2];
 		const T& operator[](uint64_t idx) const { return v[idx]; }
-		uint64_t size() const { return v.size(); }
 
-		void fill(const uint64_t size, const T num){
-			v.clear();
-			for (uint64_t i = 0; i < size; i++)
-			{
-				v.push_back(num);
-			}
+		// vector.size();
+		uint64_t size() const noexcept { return v.size(); } 
+
+
+		//###########################################
+		//#            VECTOR: == and !=            #
+		//###########################################
+		bool operator==(const Vector<T>& a) const {
+		    if (v.size() != a.size()) {
+		        return false;
+		    }
+		    for (uint64_t i = 0; i < v.size(); ++i) {
+		        if (v[i] != a[i]) {
+		            return false;
+		        }
+		    }
+		    return true;
 		}
 
-		void fill_RNG(const uint64_t size, const double min, const double max){
-			v.clear();
-			
-			std::mt19937_64 rng{}; 
-			rng.seed( std::random_device{}());
+		bool operator!=(const Vector<T>& a) const { return !(*this == a); }
 
-			for (uint64_t i = 0; i < size; i++){
+		//##################################################
+		//#            VECTOR: Scalar Addition             #
+		//##################################################
+		template <typename U, typename = typename std::enable_if<std::is_convertible<U, T>::value>::type>
+	    void inplace_add(const U a){
+	    	const uint64_t n = v.size();
+	    	for (uint64_t i = 0; i < n; ++i){
+	    		v[i] += a;
+	    	}
+	    }
+	    template <typename U, typename = typename std::enable_if<std::is_convertible<U, T>::value>::type>
+	    Vector<T> add(const U a) const{   	
+	    	Vector<T> result = *this;
+	    	result.inplace_add(a);
+	    	return result;
+	    }
+		template <typename U, typename = typename std::enable_if<std::is_convertible<U, T>::value>::type>
+		Vector<T> operator+(const U a) const {
+		    return add(a);
+		}
+		template <typename U, typename = typename std::enable_if<std::is_convertible<U, T>::value>::type>
+		Vector<T>& operator+=(const U a) {
+		    inplace_add(a);
+		    return *this;
+		}
+		//##################################################
+		//#            VECTOR: Vector Addition             #
+		//##################################################
+	    void inplace_add(const Vector<T>& a){
+	    	if (a.size() != v.size()){
+	    		throw std::runtime_error("utill:Vector.inplace_add -> Dimensions does not fit");
+	    	}
+	    	const uint64_t n = v.size();
+	    	for (uint64_t i = 0; i < n; ++i){
+	    		v[i] += a[i];
+	    	}
+	    }
+	    Vector<T> add(const Vector<T>& a) const{   	
+	    	Vector<T> result = *this;
+	    	result.inplace_add(a);
+	    	return result;
+	    }
+		Vector<T> operator+(const Vector<T>& a) const {
+		    return add(a);
+		}
+		Vector<T>& operator+=(const Vector<T>& a) {
+		    inplace_add(a);
+		    return *this;
+		}
+		//##################################################
+		//#            VECTOR: Vector Subtract             #
+		//##################################################
+	    void inplace_vec_subtract(const Vector<T>& a){
+	    	if (a.size() != v.size()){
+	    		throw std::runtime_error("utill:Vector.inplace_vec_subtract -> Dimensions does not fit");
+	    	}
 
-				v.push_back(std::uniform_real_distribution<>{min, max}(rng));
-			}
+	    	uint64_t n = a.size();
+	    	for (uint64_t i = 0; i < n; ++i){
+	    		v[i] -= a[i];
+	    	}
+
+	    }
+	    Vector<T> vec_subtract(const Vector<T>& a) const{   	
+	    	Vector<T> result = *this;
+
+	    	result.inplace_vec_subtract(a);
+
+	    	return result;
+	    }
+		Vector<T> operator-(const Vector<T>& a) const {
+		    return vec_subtract(a);
+		}
+		Vector<T>& operator-=(const Vector<T>& a) {
+		    inplace_vec_subtract(a);
+		    return *this;
+		}
+		//##################################################
+		//#            VECTOR: Vector Multiply             #
+		//##################################################
+	    void inplace_vec_multiply(const Vector<T>& a){
+	    	if (a.size() != v.size()){
+	    		throw std::runtime_error("utill:Vector.inplace_vec_multiply -> Dimensions does not fit");
+	    	}
+
+	    	uint64_t n = a.size();
+	    	for (uint64_t i = 0; i < n; ++i){
+	    		v[i] *= a[i];
+	    	}
+
+	    }
+	    Vector<T> vec_multiply(const Vector<T>& a) const{   	
+	    	Vector<T> result = *this;
+	    	result.inplace_vec_multiply(a);
+	    	return result;
+	    }
+		Vector<T> operator*(const Vector<T>& a) const {
+		    return vec_multiply(a);
+		}
+		Vector<T>& operator*=(const Vector<T>& a) {
+		    inplace_vec_multiply(a);
+		    return *this;
+		}
+		//##################################################
+		//#            VECTOR: Vector Divide               #
+		//##################################################
+	    void inplace_vec_divide(const Vector<T>& a){
+	    	if (a.size() != v.size()){
+	    		throw std::runtime_error("utill:Vector.inplace_vec_divide -> Dimensions does not fit");
+	    	}
+
+	    	uint64_t n = a.size();
+	    	for (uint64_t i = 0; i < n; ++i){
+	    		v[i] /= a[i];
+	    	}
+
+	    }
+	    Vector<T> vec_divide(const Vector<T>& a) const{   	
+	    	Vector<T> result = *this;
+	    	result.inplace_vec_divide(a);
+	    	return result;
+	    }
+		Vector<T> operator/(const Vector<T>& a) const {
+		    return vec_divide(a);
+		}
+		Vector<T>& operator/=(const Vector<T>& a) {
+		    inplace_vec_divide(a);
+		    return *this;
 		}
 
-		void disturbance(const double min, const double max){
+		//######################################################
+		//#            VECTOR: Support Functions               #
+		//######################################################
 
-			uint64_t size = v.size();
 
-			std::mt19937_64 rng{}; 
-			rng.seed(std::random_device{}());
-
-			for (uint64_t i = 0; i < size; i++){
-				v[i] *= std::uniform_real_distribution<>{min, max}(rng);
-			}
-		}
-
-		double max(){
-			double max_value = v[0];
-
-			for (uint64_t i = 1; i < v.size(); i++){
-				max_value = MAX(max_value, v[i]);
-			}
-			return max_value;
-		}
-
-		double min(){
-			double min_value = v[0];
-
-			for (uint64_t i = 1; i < v.size(); i++){
-				min_value = MIN(min_value, v[i]);
-			}
-			return min_value;
-		}
-
-		double sum() const {
-			double sum_value = 0;
-
-			for (uint64_t i = 0; i < v.size(); i++){
-				sum_value += v[i];
-			}
-			return sum_value;
-		}
-
-		double mean(){
-			double mean_value = sum()/v.size();
-			return mean_value;
-		}
 
 		inline friend std::ostream& operator << (std::ostream& out, const Vector& vec){
 			out << "[";
@@ -108,7 +213,7 @@ namespace utils{
 		void print() const{
 			std::cout << *this << std::endl;
 		}
-
+/*
 		void linspace(const T start, const T stop, const T num){
 			v.clear();
 
@@ -121,9 +226,7 @@ namespace utils{
 				}
 				//v.push_back(static_cast<T>(stop));
 			}
-			
-
-		}
+		}*/
 	};
 
 	typedef Vector<int> Vi;
@@ -133,4 +236,4 @@ namespace utils{
 
 } // namespace utils
 
-#endif // _numerics_n_
\ No newline at end of file
+#endif // _vector_n_
\ No newline at end of file
diff --git a/obj/main.o b/obj/main.o
index 64c3538..d45348f 100644
Binary files a/obj/main.o and b/obj/main.o differ
diff --git a/src/main.cpp b/src/main.cpp
index d6a6d7d..9b85861 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,45 +1,35 @@
 #include "./utils/utils.h"
 
-#include <vector>
-#include <iostream>
 
 
 
 int main(int argc, char const *argv[])
 {
 
-	utils::Vd vect1, vect2;
-	//vect1.linspace(1, 3.14159, 3);
-	vect1.fill(3, 1);
-	vect1.print();
-	//std::cout << vect[1] << std::endl;
+	utils::Vf vect1(3,1), vect2(3,1), vect3(3,1), addition_test(3,7), subtract_test(3,-2);
+	//###########################################
+	//#            VECTOR TYPE TEST             #
+	//#			   Addition Test 				#
+	//###########################################
+	vect1.inplace_add(vect2);
+	vect3 = vect1.add(vect2).add(vect1) + vect2;
+	vect3 += vect2;
+	if (vect3 != addition_test){
+		throw std::runtime_error("Addition Test -> Failed");}
 
-	
-	utils::Md matr1, matr2, matr3;
-	matr1.fill(4,4, 1);
-	matr1.print();
-	vect2 = matr1.vecmult(vect1);
-	vect2.print();
-	/*
-	matr2 = matr1.inverse("Gauss-Jordan");
-	matr3 = matr1.matmult(matr2);
-	//
-	matr1.print();
-	std::cout << "======" << std::endl;
-	matr2.print();
-	std::cout << "======" << std::endl;
-	matr3.print();
-	std::cout << "======" << std::endl;
-	matr1.matmult(matr2).print();
+	//###########################################
+	//#            VECTOR TYPE TEST             #
+	//#			   Subtract Test 				#
+	//###########################################
+	vect1.inplace_vec_subtract(vect2);
+	vect3 = vect1.vec_subtract(vect2).vec_subtract(vect1) - vect2;
+	vect3.print();
+	vect3 = vect3.add(2);
+	vect3.print();
+	if (vect3 != subtract_test){
+		throw std::runtime_error("Subtract Test -> Failed");}
 
-	std::cout << "---------" << std::endl;
-	utils::Md matr4(2,3,2);
-	vect2 = matr4.vecmult(vect1);
-	matr4.inplace_transpose();
-	matr4.transpose();
 
-	//result.print();
-	*/
 
 
 	return 0;