diff --git a/Base/Core/imstkMath.h b/Base/Core/imstkMath.h
index b64aae9432e82860a63b6620c59f3b1c9a69f923..30f1610c8a1ec4e980da33a6f3a20447ef01b576 100644
--- a/Base/Core/imstkMath.h
+++ b/Base/Core/imstkMath.h
@@ -25,6 +25,7 @@
 #include <vector>
 
 #include <Eigen/Geometry>
+#include <Eigen/Sparse>
 
 namespace imstk {
 
@@ -60,6 +61,9 @@ using Mat3d = Eigen::Matrix3d;
 using Mat4f = Eigen::Matrix4f;
 using Mat4d = Eigen::Matrix4d;
 
+// A dynamic size sparse matrix of doubles
+using SparseMatrixd = Eigen::SparseMatrix < double, Eigen::RowMajor > ;
+
 // Rigid transform (translation and rotation)
 using RigidTransform3f = Eigen::Isometry3f;
 using RigidTransform3d = Eigen::Isometry3d;
diff --git a/Base/Solvers/imstkLinearSystem.cpp b/Base/Solvers/imstkLinearSystem.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..740dc58768308d34643fe55acceb758e5e01b4de
--- /dev/null
+++ b/Base/Solvers/imstkLinearSystem.cpp
@@ -0,0 +1,126 @@
+/*=========================================================================
+
+   Library: iMSTK
+
+   Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+   & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+      http://www.apache.org/licenses/LICENSE-2.0.txt
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.
+
+   =========================================================================*/
+
+#include "imstkLinearSystem.h"
+
+#include <g3log/g3log.hpp>
+
+namespace imstk
+{
+
+template<typename SystemMatrixType>
+LinearSystem<SystemMatrixType>::LinearSystem(const SystemMatrixType &matrix, const Vectord &b)
+{
+    if (m_A.size() == m_b.size())
+    {
+        m_A = &matrix;
+        m_b = &b;
+    }
+    else
+    {
+        LOG(ERROR) << "LinearSystem::LinearSystem: The size of the matrix and the r.h.s doesn't match.";
+    }
+
+    this->F = [this](const Vectord& x) -> Vectord&
+    {
+        this->m_f = this->m_A * x;
+        return this->m_f;
+    };
+}
+
+template<typename SystemMatrixType>
+const Vectord&
+LinearSystem<SystemMatrixType>::getRHSVector() const
+{
+    return m_b;
+}
+
+template<typename SystemMatrixType>
+void
+LinearSystem<SystemMatrixType>::setRHSVector(const Vectord& newRhs)
+{
+    m_b = newRhs;
+}
+
+template<typename SystemMatrixType>
+const SystemMatrixType&
+LinearSystem<SystemMatrixType>::getMatrix() const
+{
+    return m_A;
+}
+
+template<typename SystemMatrixType>
+void
+LinearSystem<SystemMatrixType>::setMatrix(const SparseMatrixd &newMatrix)
+{
+    m_A = newMatrix;
+}
+
+template<typename SystemMatrixType>
+Vectord&
+LinearSystem<SystemMatrixType>::computeResidual(const Vectord &x, Vectord &r) const
+{
+    r = m_b - F(x);
+    return r;
+}
+
+template<typename SystemMatrixType>
+Eigen::SparseTriangularView < SystemMatrixType, Eigen::Lower >
+LinearSystem<SystemMatrixType>::getLowerTriangular() const
+{
+    return m_A.template triangularView<Eigen::Lower>();
+}
+
+template<typename SystemMatrixType>
+Eigen::SparseTriangularView < SystemMatrixType, Eigen::StrictlyLower >
+LinearSystem<SystemMatrixType>::getStrictLowerTriangular() const
+{
+    return m_A.template triangularView<Eigen::StrictlyLower>();
+}
+
+template<typename SystemMatrixType>
+Eigen::SparseTriangularView < SystemMatrixType, Eigen::Upper >
+LinearSystem<SystemMatrixType>::getUpperTrianglular() const
+{
+    return m_A.template triangularView<Eigen::Upper>();
+}
+
+template<typename SystemMatrixType>
+Eigen::SparseTriangularView < SystemMatrixType, Eigen::StrictlyUpper >
+LinearSystem<SystemMatrixType>::getStrictUpperTriangular() const
+{
+    return m_A.template triangularView<Eigen::StrictlyUpper>();
+}
+
+template<typename SystemMatrixType>
+Vectord&
+LinearSystem<SystemMatrixType>::getFunctionValue()
+{
+    return m_f;
+}
+
+template<typename SystemMatrixType> size_t
+LinearSystem<SystemMatrixType>::getSize()
+{
+    return m_A.size();
+}
+
+} //imstk
diff --git a/Base/Solvers/imstkLinearSystem.h b/Base/Solvers/imstkLinearSystem.h
new file mode 100644
index 0000000000000000000000000000000000000000..b08b6694c781d37a3dc7688eb92e321344ac30bf
--- /dev/null
+++ b/Base/Solvers/imstkLinearSystem.h
@@ -0,0 +1,119 @@
+// This file is part of the iMSTK project.
+//
+// Copyright (c) Kitware, Inc.
+//
+// Copyright (c) Center for Modeling, Simulation, and Imaging in Medicine,
+//                        Rensselaer Polytechnic Institute
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef imstkLinearSystem_h
+#define imstkLinearSystem_h
+
+#include "imstkNonlinearSystem.h"
+
+namespace imstk
+{
+
+///
+/// \class LinearSystem
+/// \brief Represents the linear system of the form \f$ Ax = b \f$
+///
+template<typename SystemMatrixType>
+class LinearSystem : public NonLinearSystem
+{
+
+public:
+    ///
+    /// \brief Constructor/destructor(s). This class can't be constructed without
+    ///     a matrix and rhs. Also, avoid copying this system.
+    ///
+    LinearSystem() = delete;
+    LinearSystem(const LinearSystem&) = delete;
+    LinearSystem& operator=(const LinearSystem&) = delete;
+    LinearSystem(const SystemMatrixType& matrix, const Vectord& b);
+
+    ///
+    /// \brief Destructor
+    ///
+    virtual ~LinearSystem() = default;
+
+    ///
+    ///  \brief Returns a reference to local right hand side vector.
+    ///
+    const Vectord& getRHSVector() const;
+
+    ///
+    /// \brief Set the system rhs corresponding to this system.
+    ///
+    void setRHSVector(const Vectord& newRhs);
+
+    ///
+    /// \brief Returns reference to local matrix.
+    ///
+    const SystemMatrixType& getMatrix() const;
+
+    ///
+    /// \brief Set the system matrix corresponding to this ODE system.
+    ///
+    void setMatrix(const SparseMatrixd& newMatrix);
+
+    ///
+    /// \brief Compute the residual as \f$\left \| b-Ax \right \|_2\f$.
+    ///
+    Vectord& computeResidual(const Vectord& x, Vectord& r) const;
+
+    ///
+    /// \brief Returns template expression for the lower triangular part of A.
+    ///
+    Eigen::SparseTriangularView < SystemMatrixType, Eigen::Lower >
+        getLowerTriangular() const;
+
+    ///
+    /// \brief Returns template expression for the strict lower triangular part of A.
+    ///
+    Eigen::SparseTriangularView < SystemMatrixType, Eigen::StrictlyLower >
+        getStrictLowerTriangular() const;
+
+    ///
+    /// \brief Returns template expression for the upper triangular part of A.
+    ///
+    Eigen::SparseTriangularView < SystemMatrixType, Eigen::Upper >
+        getUpperTrianglular() const;
+
+    ///
+    /// \brief Returns template expression for the strict upper triangular part of A.
+    ///
+    Eigen::SparseTriangularView < SystemMatrixType, Eigen::StrictlyUpper >
+        getStrictUpperTriangular() const;
+
+    ///
+    /// \brief Get the value of the function F
+    ///
+    Vectord& getFunctionValue();
+
+    ///
+    /// \brief Returns the size of the system
+    ///
+    size_t getSize();
+
+private:
+    const SystemMatrixType& m_A;
+    const Vectord& m_b;
+
+    Vectord m_f; ///> Scratch storage for matrix-vector operations
+};
+
+}
+
+#endif // imstkLinearSystem_h
diff --git a/Base/Solvers/imstkNonlinearSystem.cpp b/Base/Solvers/imstkNonlinearSystem.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6c823655b54800ddd90ab3b5a2ba77c83e81ddf
--- /dev/null
+++ b/Base/Solvers/imstkNonlinearSystem.cpp
@@ -0,0 +1,53 @@
+/*=========================================================================
+
+   Library: iMSTK
+
+   Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+   & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+      http://www.apache.org/licenses/LICENSE-2.0.txt
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.
+
+   =========================================================================*/
+
+#include "imstkNonlinearSystem.h"
+
+#include <g3log/g3log.hpp>
+
+namespace imstk
+{
+
+void
+NonLinearSystem::setFunction(const VectorFunctionType &function)
+{
+    this->m_F = function;
+}
+
+void
+NonLinearSystem::setJacobian(const MatrixFunctionType &function)
+{
+    this->m_dF = function;
+}
+
+const Vectord &
+NonLinearSystem::evaluateF(const Vectord &x)
+{
+    return this->m_F(x);
+}
+
+const SparseMatrixd &
+NonLinearSystem::evaluateJacobian(const Vectord &x)
+{
+    return this->m_dF(x);
+}
+
+} //imstk
diff --git a/Base/Solvers/imstkNonlinearSystem.h b/Base/Solvers/imstkNonlinearSystem.h
new file mode 100644
index 0000000000000000000000000000000000000000..c6055674a5b67cf8fb38b7c1e54285a747991476
--- /dev/null
+++ b/Base/Solvers/imstkNonlinearSystem.h
@@ -0,0 +1,80 @@
+// This file is part of the iMSTK project.
+//
+// Copyright (c) Kitware, Inc.
+//
+// Copyright (c) Center for Modeling, Simulation, and Imaging in Medicine,
+//                        Rensselaer Polytechnic Institute
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef imstkNonlinearSystem_h
+#define imstkNonlinearSystem_h
+
+#include <memory>
+
+#include "imstkMath.h"
+
+namespace imstk {
+
+///
+/// \class NonLinearSystem
+///
+/// \brief Base class for a nonlinear system of equations
+///
+class NonLinearSystem
+{
+public:
+    using VectorFunctionType = std::function < const Vectord&(const Vectord&) >;
+    using MatrixFunctionType = std::function < const SparseMatrixd&(const Vectord&) > ;
+
+public:
+    ///
+    /// \brief default Constructor/Destructor
+    ///
+    NonLinearSystem() = default;
+    virtual ~NonLinearSystem() = default;
+
+    ///
+    /// \brief Set function to evaluate.
+    ///
+    virtual void setFunction(const VectorFunctionType& function);
+
+    ///
+    /// \brief Set gradient function to evaluate.
+    ///
+    virtual void setJacobian(const MatrixFunctionType& function);
+
+    ///
+    /// \brief Evaluate function at specified argument
+    ///
+    /// \param x Value.
+    /// \return Function value.
+    ///
+    virtual const Vectord& evaluateF(const Vectord& x);
+
+    ///
+    /// \brief Evaluate function at specified argument
+    ///
+    /// \param x Value.
+    /// \return Function value.
+    ///
+    virtual const SparseMatrixd& evaluateJacobian(const Vectord& x);
+
+public:
+    VectorFunctionType m_F;  ///> Nonlinear function
+    MatrixFunctionType m_dF;  ///> Gradient of the Nonlinear function with respect to the unknown vector
+};
+
+}
+
+#endif // imstkNonlinearSystem_h
diff --git a/Base/TimeIntegrators/imstkTimeIntegrator.cpp b/Base/TimeIntegrators/imstkTimeIntegrator.cpp
index 6804e26dbfe80eaa1525208523039790640561eb..9f7f8604087a16bd745806d5be42c55689df2753 100644
--- a/Base/TimeIntegrators/imstkTimeIntegrator.cpp
+++ b/Base/TimeIntegrators/imstkTimeIntegrator.cpp
@@ -25,13 +25,13 @@ limitations under the License.
 
 namespace imstk {
 
-TimeIntegrator::TimeIntegrator(const Type type)
+TimeIntegrator::TimeIntegrator(const TimeIntegrator::Type type)
 {
     this->setType(type);
 }
 
 void
-TimeIntegrator::setType(const Type type)
+TimeIntegrator::setType(const TimeIntegrator::Type type)
 {
     m_type = type;
     this->setCoefficients(type);
@@ -44,19 +44,19 @@ TimeIntegrator::getType() const
 }
 
 void
-TimeIntegrator::setCoefficients(const Type type)
+TimeIntegrator::setCoefficients(const TimeIntegrator::Type type)
 {
     switch (type)
     {
-    case Type::BackwardEuler:
+    case TimeIntegrator::Type::BackwardEuler:
         m_alpha = { { 1, 0, 0 } };
         m_beta = { { 1, -1, 0 } };
         m_gamma = { { 1, -2, -1 } };
         break;
 
-    case Type::ForwardEuler:
-    case Type::NewmarkBeta:
-    case Type::CentralDifference:
+    case TimeIntegrator::Type::ForwardEuler:
+    case TimeIntegrator::Type::NewmarkBeta:
+    case TimeIntegrator::Type::CentralDifference:
         LOG(WARNING) << "TimeIntegrator::setCoefficients error: type of the time integrator not supported.";
         break;
 
diff --git a/Base/TimeIntegrators/imstkTimeIntegrator.h b/Base/TimeIntegrators/imstkTimeIntegrator.h
index 16e1a6864f7680c3dcf1b70e6dd56e84a9d071b7..1e6edf9128f28cce8fdb0c9a0b02d029909c8f25 100644
--- a/Base/TimeIntegrators/imstkTimeIntegrator.h
+++ b/Base/TimeIntegrators/imstkTimeIntegrator.h
@@ -50,7 +50,7 @@ public:
     ///
     /// \brief Constructor
     ///
-    TimeIntegrator(const Type type);
+    TimeIntegrator(const TimeIntegrator::Type type);
 
     ///
     /// \brief Destructor
@@ -60,16 +60,16 @@ public:
     ///
     /// \brief Set/Get type of the time integrator
     ///
-    void setType(const Type type);
+    void setType(const TimeIntegrator::Type type);
     const Type& getType() const;
 
     ///
     /// \brief Set coefficients for a given time integrator type
     ///
-    void setCoefficients(const Type type);
+    void setCoefficients(const TimeIntegrator::Type type);
 
 protected:
-    Type m_type; ///> Type of the time integrator
+    TimeIntegrator::Type m_type; ///> Type of the time integrator
 
     // Coefficients of the time integrator
     std::array<double,3> m_alpha;
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8b77ef8f3a80dc44ea7a9dcb16735e9572c1eb58..299fad1f0bcff05483c56dab374bfac5076c9f15 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -161,6 +161,8 @@ add_subdirectory(Base/Collision)
 add_subdirectory(Base/Scene)
 add_subdirectory(Base/Rendering)
 add_subdirectory(Base/SimulationManager)
+add_subdirectory(Base/Solvers)
+add_subdirectory(Base/TimeIntegrators)
 
 #--------------------------------------------------------------------------
 # Export Targets