diff --git a/bin/abc_lab b/bin/abc_lab
new file mode 100755
index 0000000..ac1a151
Binary files /dev/null and b/bin/abc_lab differ
diff --git a/include/core/global_config.h b/include/core/global_config.h
new file mode 100644
index 0000000..c18e981
--- /dev/null
+++ b/include/core/global_config.h
@@ -0,0 +1,36 @@
+#pragma once
+
+
+namespace core{
+
+    enum class GridKind { Uniform, NonUniform };
+    enum class FDKind   { Central, Forward, Backward };
+    enum class BCKind   { Dirichlet, Neumann /*, Robin*/ };
+    enum class SolverKind { LU, Inverse /*, CG*/ };
+
+    template <typename T>
+    struct BC { 
+        FDKind fd{FDKind::Forward};
+        BCKind kind{BCKind::Dirichlet}; 
+        T value{T(0)}; 
+    };
+
+
+    // Global default config holder
+    template <typename T>
+    struct Configs {
+        GridKind grid{GridKind::Uniform};
+        FDKind   fd{FDKind::Central};
+        BC<T>    left{FDKind::Forward, BCKind::Dirichlet, T(0) };
+        BC<T>    right{FDKind::Backward, BCKind::Dirichlet, T(0) };
+        SolverKind solver{SolverKind::LU};
+
+        static Configs& defaults() {
+            static Configs g{};   // process-wide defaults
+            return g;
+        }
+    };
+
+
+
+} // namespace core
diff --git a/include/core/omp_config.h b/include/core/omp_config.h
index 91b7e75..6449035 100644
--- a/include/core/omp_config.h
+++ b/include/core/omp_config.h
@@ -1,6 +1,6 @@
-
-
 #pragma once
+
+#include <vector>
 #include <omp.h>
 
 
diff --git a/include/modules/fluids/diffusion1d.h b/include/modules/fluids/diffusion1d.h
index e69de29..221af99 100644
--- a/include/modules/fluids/diffusion1d.h
+++ b/include/modules/fluids/diffusion1d.h
@@ -0,0 +1,129 @@
+#pragma once
+
+#include "modules/mesh/mesh1d.h"
+#include "core/global_config.h"
+#include "utils/matrix.h"
+#include "utils/vector.h"
+
+
+
+
+namespace fluids {
+
+	template <typename T>
+	struct Diffusion1D{
+		const core::Configs<T>& cfg;
+		const mesh::Mesh1D<T>& mesh;
+		T Gamma{1};
+
+	
+
+		// Constructor
+		Diffusion1D(const core::Configs<T>& configs, const mesh::Mesh1D<T>& Mesh, T Gamma_const=T(1)): cfg(configs), mesh(Mesh), Gamma(Gamma_const) {}
+
+
+
+
+
+
+		void assemble(utils::Matrix<T>& A, utils::Vector<T>& b, utils::Vector<T>& s){
+
+			uint64_t N = mesh.center_idx + 1;
+
+			if (N < 3){
+				throw std::runtime_error("Diffusion1D: need N>=3");
+			}
+			if (A.rows() != N || A.cols() != N){
+				A = utils::Matrix<T>(N, N, T(0));
+			}
+			if (b.size() != N){
+				b = utils::Vector<T>(N, T(0));
+			}
+		
+
+			if (cfg.grid == core::GridKind::Uniform){
+				// Core of A
+				if (cfg.fd == core::FDKind::Central){
+					uniform_central_finite_diffrence_2_order(A,b,s);
+				}
+
+
+
+				// Left BC of A
+				if (cfg.left.kind == core::BCKind::Dirichlet){
+					BC_uniform_backward_finite_diffrence_2_order_Dirichlet(A,b,s);
+				}else if (cfg.left.kind == core::BCKind::Neumann){
+					BC_uniform_backward_finite_diffrence_2_order_Neumann(A,b,s);
+				}
+			}
+		}
+
+
+		void uniform_central_finite_diffrence_2_order(utils::Matrix<T>& A, utils::Vector<T>& b, utils::Vector<T>& s){
+			T xm, xc, xp;
+
+			for (uint64_t i = 1; i < mesh.center_idx; ++i){
+				xm = mesh.center(i-1);
+				xc = mesh.center(i);
+				xp = mesh.center(i+1);
+
+				A(i, i-1) = Gamma/(xc - xm);
+				A(i, i) = -((Gamma/(xp - xc)) + (Gamma/(xc - xm)));
+				A(i, i+1) = Gamma/(xp - xc);
+
+				b[i] = -s[i]*mesh.dx(i);
+			}
+		}
+
+		void BC_uniform_backward_finite_diffrence_2_order_Dirichlet(utils::Matrix<T>& A, utils::Vector<T>& b, utils::Vector<T>& s){
+			T xm;
+			T xw = mesh.vertice(0);
+			T xc = mesh.center(0);
+			T xp = mesh.center(1);
+			T xe;
+			uint64_t N = mesh.center_idx;
+
+			A(0, 0) = -((Gamma/(xp - xc)) + (Gamma/(xc - xw)));
+			A(0, 1) = Gamma/(xp - xc);
+
+			b[0] = -s[0]*mesh.dx(0) - Gamma*(cfg.left.value/(xc - xw));
+
+			xm = mesh.center(N-1);
+			xc = mesh.center(N);
+			xe = mesh.vertice(N+1);
+
+			A(N, N-1) = Gamma/(xc - xm);
+			A(N, N) = -((Gamma/(xe - xc)) + (Gamma/(xc - xm)));
+
+			b[N] = -s[N]*mesh.dx(N) - Gamma*(cfg.right.value/(xe - xc));
+			A.print();
+			b.print();
+		}
+
+		void BC_uniform_backward_finite_diffrence_2_order_Neumann(utils::Matrix<T>& A, utils::Vector<T>& b, utils::Vector<T>& s){
+			T xm;
+			T xc = mesh.center(0);
+			T xp = mesh.center(1);
+
+			uint64_t N = mesh.center_idx;
+
+			A(0, 0) = -Gamma/(xp - xc);
+			A(0, 1) = Gamma/(xp - xc);
+
+			b[0] = -s[0]*mesh.dx(0) - (Gamma*cfg.left.value);
+
+			xm = mesh.center(N-1);
+			xc = mesh.center(N);
+
+			A(N, N-1) = Gamma/(xc - xm);
+			A(N, N) = -Gamma/(xc - xm);
+
+			b[N] = -s[N]*mesh.dx(N) - Gamma*(cfg.right.value);
+
+			A.print();
+			b.print();
+		}
+
+	};
+
+} // end namespace fluids
\ No newline at end of file
diff --git a/include/modules/fluids/fluids.h b/include/modules/fluids/fluids.h
new file mode 100644
index 0000000..066af5f
--- /dev/null
+++ b/include/modules/fluids/fluids.h
@@ -0,0 +1,5 @@
+#pragma once
+
+#include "modules/fluids/diffusion1d.h"
+
+
diff --git a/include/modules/mesh/mesh.h b/include/modules/mesh/mesh.h
new file mode 100644
index 0000000..369ebd7
--- /dev/null
+++ b/include/modules/mesh/mesh.h
@@ -0,0 +1,4 @@
+#pragma once
+
+#include "modules/mesh/mesh1d.h"
+
diff --git a/include/modules/grid1d.h b/include/modules/mesh/mesh1d.h
similarity index 77%
rename from include/modules/grid1d.h
rename to include/modules/mesh/mesh1d.h
index 8fc35fc..5373580 100644
--- a/include/modules/grid1d.h
+++ b/include/modules/mesh/mesh1d.h
@@ -2,20 +2,18 @@
 
 #include "utils/vector.h"
 
-
-namespace fvm {
+namespace mesh {
 
 	template <typename T>
-
-	struct Grid1D{
+	struct Mesh1D{
 		uint64_t center_idx;				// max cell index
 		uint64_t vertices_idx;				// max vertice index
 		utils::Vector<T> centers;			// size N (unknowns at cell centers)
 		utils::Vector<T> vertices;			// size N+1 (face positions)
 
-		Grid1D() = default;
+		Mesh1D() = default;
 
-		explicit Grid1D(const utils::Vector<T>& midpoints){
+		explicit Mesh1D(const utils::Vector<T>& midpoints){
 			centers = midpoints;
 			center_idx = centers.size()-1;
 			vertices_idx = centers.size();		
@@ -33,11 +31,12 @@ namespace fvm {
 
 		T dx(uint64_t i) const { check(i); return vertices[i+1] - vertices[i]; }
 		T center(uint64_t i) const { check(i); return centers[i]; }
+		T vertice(uint64_t i) const {; return vertices[i]; }
 
 	    void check(uint64_t i) const {
-	        if (i > center_idx) throw std::runtime_error("Grid1D: cell index out of range");
+	        if (i > center_idx) throw std::runtime_error("Mesh1D: cell index out of range");
 	    }
 	};
 
 
-}
\ No newline at end of file
+} // end namespace mesh
\ No newline at end of file
diff --git a/include/numerics/numerics.h b/include/numerics/numerics.h
index e3cd5ad..bb4760e 100644
--- a/include/numerics/numerics.h
+++ b/include/numerics/numerics.h
@@ -10,10 +10,6 @@
 #include "./numerics/min.h"
 #include "./numerics/max.h"
 #include "./numerics/abs.h"
-#include "./numerics/interpolation1d_base.h"                       		// base
-#include "./numerics/interpolation1d/interpolation1d_linear.h"  		// derived
-#include "./numerics/interpolation1d/interpolation1d_polynomial.h"  	// derived
-#include "./numerics/interpolation1d/interpolation1d_cubic_spline.h"  	// derived
-#include "./numerics/interpolation1d/interpolation1d_rational.h"  		// derived 
-#include "./numerics/interpolation1d/interpolation1d_barycentric.h"  		// derived
+#include "./numerics/interpolation1d.h"                       		// base
+
 
diff --git a/include/utils/generators.h b/include/utils/generators.h
new file mode 100644
index 0000000..3369d29
--- /dev/null
+++ b/include/utils/generators.h
@@ -0,0 +1,4 @@
+// "./utils/generators.h"
+#pragma once
+
+#include "./utils/generators/linspace.h"
diff --git a/include/utils/generators/linspace.h b/include/utils/generators/linspace.h
new file mode 100644
index 0000000..9eb5086
--- /dev/null
+++ b/include/utils/generators/linspace.h
@@ -0,0 +1,37 @@
+#pragma once
+
+#include "utils/vector.h"
+
+
+namespace utils{
+
+	template <typename T>
+	void inplace_linspace(utils::Vector<T>& a, T start, T stop, bool endpoint=true){
+
+		uint64_t N = a.size();
+		T step;
+
+		if (endpoint){
+			step = (stop - start) / static_cast<T>(N - 1);
+		}else{
+			step = (stop - start) / static_cast<T>(N);
+		}
+
+		for (uint64_t i = 0; i < N; ++i){
+			a[i] = start + (step*static_cast<T>(i));
+		}
+	}
+
+	template <typename T>
+	utils::Vector<T> linspace(T start, T stop, uint64_t N, bool endpoint=true){
+
+		utils::Vector<T> a(N);
+
+		inplace_linspace(a, start, stop, endpoint);
+
+		return a;
+	}
+
+
+
+} // end namespace utils
\ No newline at end of file
diff --git a/include/utils/utils.h b/include/utils/utils.h
index cfd11c5..86904a4 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/generators.h"
\ No newline at end of file
diff --git a/obj/main.d b/obj/main.d
new file mode 100644
index 0000000..be73496
--- /dev/null
+++ b/obj/main.d
@@ -0,0 +1,56 @@
+obj/main.o: src/main.cpp include/./core/omp_config.h \
+ include/./utils/utils.h include/./utils/vector.h \
+ include/./utils/matrix.h include/./utils/generators.h \
+ include/./utils/generators/linspace.h include/utils/vector.h \
+ include/./numerics/numerics.h include/./numerics/initializers/eye.h \
+ include/./numerics/matequal.h include/./numerics/abs.h \
+ include/./numerics/transpose.h include/./numerics/inverse.h \
+ include/./numerics/inverse/inverse_gauss_jordan.h \
+ include/./numerics/inverse/inverse_lu.h include/./decomp/lu.h \
+ include/./numerics/matmul.h include/./numerics/matvec.h \
+ include/./numerics/min.h include/./numerics/max.h \
+ include/./numerics/interpolation1d.h \
+ include/./numerics/interpolation1d/interpolation1d_barycentric.h \
+ include/./numerics/interpolation1d/interpolation1d_base.h \
+ include/./numerics/interpolation1d/interpolation1d_cubic_spline.h \
+ include/./numerics/interpolation1d/interpolation1d_linear.h \
+ include/./numerics/interpolation1d/interpolation1d_polynomial.h \
+ include/./numerics/interpolation1d/interpolation1d_rational.h \
+ include/./decomp/decomp.h include/./modules/mesh/mesh.h \
+ include/modules/mesh/mesh1d.h include/modules/fluids/fluids.h \
+ include/modules/fluids/diffusion1d.h include/core/global_config.h \
+ include/utils/matrix.h
+include/./core/omp_config.h:
+include/./utils/utils.h:
+include/./utils/vector.h:
+include/./utils/matrix.h:
+include/./utils/generators.h:
+include/./utils/generators/linspace.h:
+include/utils/vector.h:
+include/./numerics/numerics.h:
+include/./numerics/initializers/eye.h:
+include/./numerics/matequal.h:
+include/./numerics/abs.h:
+include/./numerics/transpose.h:
+include/./numerics/inverse.h:
+include/./numerics/inverse/inverse_gauss_jordan.h:
+include/./numerics/inverse/inverse_lu.h:
+include/./decomp/lu.h:
+include/./numerics/matmul.h:
+include/./numerics/matvec.h:
+include/./numerics/min.h:
+include/./numerics/max.h:
+include/./numerics/interpolation1d.h:
+include/./numerics/interpolation1d/interpolation1d_barycentric.h:
+include/./numerics/interpolation1d/interpolation1d_base.h:
+include/./numerics/interpolation1d/interpolation1d_cubic_spline.h:
+include/./numerics/interpolation1d/interpolation1d_linear.h:
+include/./numerics/interpolation1d/interpolation1d_polynomial.h:
+include/./numerics/interpolation1d/interpolation1d_rational.h:
+include/./decomp/decomp.h:
+include/./modules/mesh/mesh.h:
+include/modules/mesh/mesh1d.h:
+include/modules/fluids/fluids.h:
+include/modules/fluids/diffusion1d.h:
+include/core/global_config.h:
+include/utils/matrix.h:
diff --git a/obj/main.o b/obj/main.o
new file mode 100644
index 0000000..15bd227
Binary files /dev/null and b/obj/main.o differ
diff --git a/src/main.cpp b/src/main.cpp
index 0d20d31..5efbf29 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,146 +1,46 @@
+#include "./core/omp_config.h"
+
 #include "./utils/utils.h"
 #include "./numerics/numerics.h"
 #include "./decomp/decomp.h"
-#include "./core/omp_config.h"
-
-#include "./modules/grid1d.h"
-#include "./modules/finitedifference1d.h"
 
 
-#include <iostream>
-#include <stdexcept>
+#include "./modules/mesh/mesh.h"
+#include "modules/fluids/fluids.h"
 
 
-#include <chrono>
 
-//#include "./numerics/interpolation/interpolation_linear.h"
+//#include <iostream>
+//#include <stdexcept>
+//#include <chrono>
+
 
 
 
 
 int main(int argc, char const *argv[])
 {    
-    utils::Md A;
+    utils::Vd a = utils::linspace<double>(1, 10, 10, true);
+    a.print();
+    mesh::Mesh1D<double> mesh(a);
+    mesh.generate_vertices(0.5, 10.5);
+    double Gamma = 1.0;
+
     
-/*
-    int hw = omp_get_max_active_levels();
-    if (hw <= 1) return 0;
-
-    const uint64_t m=512, k=512, p=512; // ~134M MACs; adjust if needed
-    utils::Md A(m,k,1), B(k,p,1), C(k,p,1);
-
-    omp_set_max_active_levels(1);
-
-    auto t0 = std::chrono::high_resolution_clock::now();
-    for (int i = 0; i < m*k*p; ++i){
-        A==B
-    }
-    double t1 = std::chrono::duration<double>(std::chrono::high_resolution_clock::now() - t0).count();
+    utils::Md A;
+    utils::Vd b, s(10,1);
 
 
-    omp_set_max_active_levels(2);
-    auto t0 = std::chrono::high_resolution_clock::now();
-    for (int i = 0; i < m*k*p; ++i){
-        A==B
-    }
-    double t1 = std::chrono::duration<double>(std::chrono::high_resolution_clock::now() - t0).count();
-
-    omp_set_num_threads(prev);
-
-    // Must not be notably slower with many threads
-    CHECK(tN <= t1 * 1.05, "rows_omp: multi-thread slower than single-thread");
+    core::Configs<double>& cfg = core::Configs<double>::defaults();
+    cfg.grid = core::GridKind::Uniform;
+    cfg.left = {core::FDKind::Forward, core::BCKind::Neumann, 0.0};
+    cfg.right = {core::FDKind::Backward, core::BCKind::Neumann, 0.0};
+    cfg.solver    = core::SolverKind::LU;
 
 
+    fluids::Diffusion1D<double> diffusion(cfg, mesh, Gamma);
+    diffusion.assemble(A, b, s);
 
-    utils::Md A(5,5, 1);
-    utils::Md B(5,5, 1);
-    utils::Md C(5,5, 2);
-
-    bool result1 = (A==B);
-    bool result2 = (A==C);
-
-    omp_set_max_active_levels(1):
-
-    for (int i = 0; i < 100; ++i){
-        (A==B)
-    }
-
-    omp_set_max_active_levels(2):
-
-    for (int i = 0; i < 100; ++i){
-        (A==B)
-    }
-
-    std::cout << result1 << std::endl;
-    std::cout << result2 << std::endl;
-
-
-
-
-*/
-
-
-    /*
-
-    utils::Vector<double> x(100, 0), y(100,0);
-    for (uint64_t i = 0; i < 100; ++i){
-        x[i] = i+1;
-        y[i] = 1 + i*0.1;
-    }
-    //y[9] = 1.5;
-
-    x.print();
-    y.print();
-
-    double p = 5.5;
-
-
-
-    numerics::interp_linear<double> linear(x,y);
-    numerics::interp_polynomial<double> polynomial(x,y, 3);
-    numerics::interp_cubic_spline<double> cubic_spline(x,y);
-    numerics::interp_rational<double> rational(x,y,2);
-    numerics::interp_barycentric<double> barycentric(x,y, 2);
-
-    std::cout << "=== interpolate: " << p << " ===" << std::endl;
-
-    std::cout << linear.interp(p) << std::endl;
-    std::cout << linear.interp(p) << std::endl;
-    std::cout << polynomial.interp(p) << std::endl; // error = polynomial.dy
-    std::cout << cubic_spline.interp(p) << std::endl;
-    std::cout << rational.interp(p) << std::endl;
-    std::cout << barycentric.interp(p) << std::endl;
-
-    p += 0.01;
-
-    std::cout << "=== interpolate: " << p << " ===" << std::endl;
-
-    std::cout << linear.interp(p) << std::endl;
-    std::cout << polynomial.interp(p) << std::endl;
-    std::cout << cubic_spline.interp(p) << std::endl;
-    std::cout << rational.interp(p) << std::endl;
-    std::cout << barycentric.interp(p) << std::endl;
-
-    p += 0.01;
-
-    std::cout << "=== interpolate: " << p << " ===" << std::endl;
-
-    std::cout << linear.interp(p) << std::endl;
-    std::cout << polynomial.interp(p) << std::endl;
-    std::cout << cubic_spline.interp(p) << std::endl;
-    std::cout << rational.interp(p) << std::endl;
-    std::cout << barycentric.interp(p) << std::endl;
-
-    p += 50.01;
-
-    std::cout << "=== interpolate: " << p << " ===" << std::endl;
-
-    std::cout << linear.interp(p) << std::endl;
-    std::cout << polynomial.interp(p) << std::endl;
-    std::cout << cubic_spline.interp(p) << std::endl;
-    std::cout << rational.interp(p) << std::endl;
-    std::cout << barycentric.interp(p) << std::endl;
-
-*/
+    
     return 0;
 }
\ No newline at end of file