diff --git a/CHANGELOG.md b/CHANGELOG.md
index a76d2e095a..063f32b16a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -48,6 +48,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Improved `xmldocchecker` with XSD schema validation for XML files.
 
+- [#1458](http://github.com/nelson-lang/nelson/issues/1458): Eigen 5.0.0 used if detected.
+
 ### Fixed
 
 - `ans` variable: Fixed issue to ensure 'ans' is only created for expressions, preventing unintended variable creation in the workspace.
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0b1a2d6c40..39a7b1880b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -244,7 +244,20 @@ if(${CMAKE_SYSTEM_NAME} MATCHES "Linux")
   endif()
 endif()
 # ==============================================================================
-find_package(Eigen3 3.3 REQUIRED)
+find_package(Eigen3 3.3 QUIET NO_MODULE)
+if (NOT Eigen3_FOUND)
+  find_package(Eigen3 REQUIRED)
+endif()
+if (TARGET Eigen3::Eigen)
+  message(STATUS "Found Eigen3 target: Eigen3::Eigen")
+  get_target_property(EIGEN3_INCLUDE_DIR Eigen3::Eigen INTERFACE_INCLUDE_DIRECTORIES)
+else()
+  # In case we found Eigen3 but not as a target
+  if (NOT DEFINED EIGEN3_INCLUDE_DIR AND DEFINED EIGEN3_INCLUDE_DIRS)
+    set(EIGEN3_INCLUDE_DIR ${EIGEN3_INCLUDE_DIRS})
+  endif()
+endif()
+message(STATUS "Eigen3 include directory: ${EIGEN3_INCLUDE_DIR}")
 # ==============================================================================
 if(NOT WITHOUT_GUI_MODULE)
   include(QtConfig)
diff --git a/modules/operators/src/cpp/LeftDivide.cpp b/modules/operators/src/cpp/LeftDivide.cpp
index 7ef7c4ec1b..5c340c79d6 100644
--- a/modules/operators/src/cpp/LeftDivide.cpp
+++ b/modules/operators/src/cpp/LeftDivide.cpp
@@ -19,41 +19,6 @@
 //=============================================================================
 namespace Nelson {
 //=============================================================================
-// static bool
-// promoteCommonType(ArrayOf& A, ArrayOf& B)
-//{
-//    bool wasPromoted = true;
-//    if (A.getDataClass() != B.getDataClass()) {
-//        if (A.isComplex() || B.isComplex()) {
-//            if (A.isComplex()) {
-//                if (B.getDataClass() == NLS_DOUBLE) {
-//                    B.promoteType(NLS_DCOMPLEX);
-//                } else if (B.getDataClass() == NLS_SINGLE) {
-//                    B.promoteType(NLS_SCOMPLEX);
-//                } else {
-//                    wasPromoted = false;
-//                }
-//            } else {
-//                if (A.getDataClass() == NLS_DOUBLE) {
-//                    A.promoteType(NLS_DCOMPLEX);
-//                } else if (A.getDataClass() == NLS_SINGLE) {
-//                    A.promoteType(NLS_SCOMPLEX);
-//                } else {
-//                    wasPromoted = false;
-//                }
-//            }
-//        } else {
-//            if (A.getDataClass() == NLS_SINGLE || B.getDataClass() == NLS_SINGLE) {
-//                A.promoteType(NLS_SINGLE);
-//                B.promoteType(NLS_SINGLE);
-//            } else {
-//                wasPromoted = false;
-//            }
-//        }
-//    }
-//    return wasPromoted;
-//}
-////=============================================================================
 ArrayOf
 LeftDivide(ArrayOf A, ArrayOf B, bool& needToOverload)
 {
diff --git a/modules/operators/src/cpp/SVDDecompositionComplexTemplate.hpp b/modules/operators/src/cpp/SVDDecompositionComplexTemplate.hpp
index 2d5c1ffeaf..28549452b0 100644
--- a/modules/operators/src/cpp/SVDDecompositionComplexTemplate.hpp
+++ b/modules/operators/src/cpp/SVDDecompositionComplexTemplate.hpp
@@ -38,10 +38,9 @@ solveSVDDecompositionComplex(NelsonType destinationClass, const ArrayOf& matA, c
     Eigen::Map<Eigen::Matrix<std::complex<T>, Eigen::Dynamic, Eigen::Dynamic>> matCz(
         ptrCz, outDim.getRows(), outDim.getColumns());
 
-    matCz = matAz
-                .bdcSvd(Eigen::DecompositionOptions::ComputeThinU
-                    | Eigen::DecompositionOptions::ComputeThinV)
-                .solve(matBz);
+    Eigen::JacobiSVD<Eigen::Matrix<std::complex<T>, Eigen::Dynamic, Eigen::Dynamic>> svd(
+        matAz, Eigen::ComputeThinU | Eigen::ComputeThinV);
+    matCz = svd.solve(matBz);
     return res;
 }
 //=============================================================================
diff --git a/modules/operators/src/cpp/SVDDecompositionRealTemplate.hpp b/modules/operators/src/cpp/SVDDecompositionRealTemplate.hpp
index 2f1cebaac9..0bfdfd8402 100644
--- a/modules/operators/src/cpp/SVDDecompositionRealTemplate.hpp
+++ b/modules/operators/src/cpp/SVDDecompositionRealTemplate.hpp
@@ -34,10 +34,9 @@ solveSVDDecompositionReal(NelsonType destinationClass, const ArrayOf& matA, cons
     Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>> matCz(
         ptrC, outDim.getRows(), outDim.getColumns());
 
-    matCz = matAz
-                .bdcSvd(Eigen::DecompositionOptions::ComputeThinU
-                    | Eigen::DecompositionOptions::ComputeThinV)
-                .solve(matBz);
+    Eigen::JacobiSVD<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>> svd(
+        matAz, Eigen::ComputeFullU | Eigen::ComputeFullV);
+    matCz = svd.solve(matBz);
     return res;
 }
 //=============================================================================