diff --git a/ALFI/ALFI/spline/cubic.h b/ALFI/ALFI/spline/cubic.h
index c290f23..cd31df8 100644
--- a/ALFI/ALFI/spline/cubic.h
+++ b/ALFI/ALFI/spline/cubic.h
@@ -303,7 +303,7 @@ namespace alfi::spline {
 						const auto h12 = X[2]-X[0], h13 = X[3]-X[0];
 						const auto d1 = Y[1]-Y[0], d12 = Y[2]-Y[0], d13 = Y[3]-Y[0];
 
-						const auto abc = util::linalg::linsolve(
+						const auto abc = util::linalg::lup_solve(
 							{{std::pow(h1, 3), std::pow(h1, 2), h1},
 								{std::pow(h12, 3), std::pow(h12, 2), h12},
 								{std::pow(h13, 3), std::pow(h13, 2), h13}},
@@ -350,7 +350,7 @@ namespace alfi::spline {
 						const auto h12 = X[n-2]-X[n-4], h13 = X[n-1]-X[n-4];
 						const auto d1 = dY[n-4], d12 = Y[n-2] - Y[n-4], d13 = Y[n-1] - Y[n-4];
 
-						const auto abc = util::linalg::linsolve(
+						const auto abc = util::linalg::lup_solve(
 							{{std::pow(h1, 3), std::pow(h1, 2), h1},
 								{std::pow(h12, 3), std::pow(h12, 2), h12},
 								{std::pow(h13, 3), std::pow(h13, 2), h13}},
diff --git a/ALFI/ALFI/util/linalg.h b/ALFI/ALFI/util/linalg.h
index 95410c7..1cf23b2 100644
--- a/ALFI/ALFI/util/linalg.h
+++ b/ALFI/ALFI/util/linalg.h
@@ -5,62 +5,73 @@
 #include "../config.h"
 
 namespace alfi::util::linalg {
-	template <typename Number = DefaultNumber, template <typename> class Container = DefaultContainer>
-	Container<Number> linsolve(Container<Container<Number>>&& A, const Container<Number>& B, Number epsilon = std::numeric_limits<Number>::epsilon()) {
-		const auto N = B.size();
+	/**
+		@brief Solves a system of linear equations using LUP decomposition.
 
-		Container<SizeT> P(N);
-		for (SizeT i = 0; i < N; ++i) {
-			P[i] = i;
-		}
+		This function solves the system of linear equations \f(AX = B\f), where \f(A\f) is a square matrix,
+		by performing LUP decomposition and forward-backward substitution.\n
+		The input matrix \f(A\f) is decomposed into lower \f(L\f) and upper \f(U\f) triangular matrices,
+		such that \f(PA = LU\f), transforming the equation into \f(PLUX = B\f) or \f(LUX = PB\f).
+
+		The solution is found in two stages:
+		1. The first stage solves the system \f(LY = PB\f) for the intermediate vector \f(Y\f) using forward substitution.
+		2. The second stage solves the system \f(UX = Y\f) for the final solution \f(X\f) using backward substitution.
+
+		If the pivot element in a given column is too small (less than the specified \f(epsilon\f)),
+		the function returns an empty container, indicating that the matrix is degenerate,
+		and prints the corresponding message to `std::cerr`.
 
-		for (SizeT i = 0; i < N; ++i) {
+		The input matrix \f(A\f) and vector \f(B\f) are modified in place.
+
+		@param A the matrix (2D container of size \f(n \times n\f))
+		@param B the right-hand side vector (1D container of size \f(n\f))
+		@param epsilon a small threshold value to detect degeneracy (default is machine epsilon)
+		@return the solution vector or an empty container if the matrix is degenerate
+	 */
+	template <typename Number = DefaultNumber, template <typename> class Container = DefaultContainer>
+	Container<Number> lup_solve(Container<Container<Number>>&& A, Container<Number>&& B, Number epsilon = std::numeric_limits<Number>::epsilon()) {
+		const auto n = B.size();
+		assert(n == A.size());
+		for (SizeT i = 0; i < n; ++i) {
+			assert(n == A[i].size());
 			Number max_a = std::abs(A[i][i]);
 			SizeT i_max = i;
-
-			for (SizeT k = i + 1; k < N; ++k) {
+			for (SizeT k = i + 1; k < n; ++k) {
 				const auto cur = std::abs(A[k][i]);
 				if (cur > max_a) {
 					max_a = cur;
 					i_max = k;
 				}
 			}
-
 			if (max_a < epsilon) {
 				std::cerr << "Error in function " << __FUNCTION__ << ": Matrix A is degenerate. Returning an empty array..." << std::endl;
 				return {};
 			}
-
 			if (i_max != i) {
-				std::swap(P[i], P[i_max]);
 				std::swap(A[i], A[i_max]);
+				std::swap(B[i], B[i_max]);
 			}
-
-			for (SizeT j = i + 1; j < N; ++j) {
+			for (SizeT j = i + 1; j < n; ++j) {
 				A[j][i] /= A[i][i];
-				for (SizeT k = i + 1; k < N; ++k) {
+				for (SizeT k = i + 1; k < n; ++k) {
 					A[j][k] -= A[j][i] * A[i][k];
 				}
 			}
 		}
-
-		Container<Number> X(N);
-
-		for (SizeT i = 0; i < N; ++i) {
-			X[i] = B[P[i]];
+		Container<Number> X(n);
+		for (SizeT i = 0; i < n; ++i) {
+			X[i] = B[i];
 			for (SizeT k = 0; k < i; ++k) {
 				X[i] -= A[i][k] * X[k];
 			}
 		}
-
-		for (SizeT iter = 0; iter < N; ++iter) {
-			const auto i = N - 1 - iter;
-			for (SizeT k = i + 1; k < N; ++k) {
+		for (SizeT iter = 0; iter < n; ++iter) {
+			const auto i = n - 1 - iter;
+			for (SizeT k = i + 1; k < n; ++k) {
 				X[i] -= A[i][k] * X[k];
 			}
 			X[i] /= A[i][i];
 		}
-
 		return X;
 	}
 
diff --git a/ALFI/ALFI/util/spline.h b/ALFI/ALFI/util/spline.h
index 49396ff..c9ab1ad 100644
--- a/ALFI/ALFI/util/spline.h
+++ b/ALFI/ALFI/util/spline.h
@@ -65,7 +65,7 @@ namespace alfi::util::spline {
 			const auto h12 = X[2]-X[0], h13 = X[3]-X[0];
 			const auto d1 = Y[1]-Y[0], d12 = Y[2]-Y[0], d13 = Y[3]-Y[0];
 
-			const auto abc = util::linalg::linsolve(
+			const auto abc = util::linalg::lup_solve(
 				{{std::pow(h1, 3), std::pow(h1, 2), h1},
 					{std::pow(h12, 3), std::pow(h12, 2), h12},
 					{std::pow(h13, 3), std::pow(h13, 2), h13}},