diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 7e89028aa145..4f8bd2079fd1 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -1222,7 +1222,8 @@ class CConfig {
   Supercatalytic_Wall;                      /*!< \brief Flag for supercatalytic wall. */
   string GasModel,                          /*!< \brief Gas Model. */
   *Wall_Catalytic;                          /*!< \brief Pointer to catalytic walls. */
-  string NoneqStateModel;                   /*!< \brief Nonequilibtrium State Model. */
+  string NoneqStateModel;                   /*!< \brief Nonequilibrium State Model. */
+  bool Newton_2T;                           /*!< \brief M++ derivative vs. Perturbation Method for Newton's Method. */
   TRANSCOEFFMODEL   Kind_TransCoeffModel;   /*!< \brief Transport coefficient Model for NEMO solver. */
   su2double CatalyticEfficiency;            /*!< \brief Wall catalytic efficiency. */
   su2double *Inlet_MassFrac;                /*!< \brief Specified Mass fraction vectors for NEMO inlet boundaries. */
@@ -3822,6 +3823,13 @@ class CConfig {
    */
   string GetNoneqStateModel(void) const {return NoneqStateModel;}
 
+    /*!
+   * \brief 2T Nonequilibrium numerical method to solve for temperatures.
+   * \return Bool to use M++ implementation vs. Perturbation Method.
+   */
+  bool Get2TNewton(void) const {return Newton_2T;}
+
+
   /*!
    * \brief Get the transport coefficient model.
    * \return Index of transport coefficient model.
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index a3609eabc4ea..6de8b97a1b9f 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -1228,8 +1228,10 @@ void CConfig::SetConfig_Options() {
   /*--- Reading gas model as string or integer depending on TC library used. ---*/
   /* DESCRIPTION: Specify chemical model for multi-species simulations - read by Mutation++ library*/
   addStringOption("GAS_MODEL", GasModel, string("N2"));
-  /* DESCRIPTION: Specify transport coefficient model for multi-species simulations */
+  /* DESCRIPTION: Specify Nonequilibrium Temperature Model */
   addStringOption("NONEQ_STATE_MODEL", NoneqStateModel, string("2T"));
+  /* DESCRIPTION: Specify M++ Derivative vs. Perturbation Method for Newton's Method solving for temperatures */
+  addBoolOption("2T_NEWTON", Newton_2T, true);
   /* DESCRIPTION: Specify transport coefficient model for multi-species simulations */
   addEnumOption("TRANSPORT_COEFF_MODEL", Kind_TransCoeffModel, TransCoeffModel_Map, TRANSCOEFFMODEL::WILKE);
   /* DESCRIPTION: Specify mass fraction of each species */
diff --git a/SU2_CFD/include/fluid/CMutationTCLib.hpp b/SU2_CFD/include/fluid/CMutationTCLib.hpp
index f1799579b4a1..7a25ca1eda85 100644
--- a/SU2_CFD/include/fluid/CMutationTCLib.hpp
+++ b/SU2_CFD/include/fluid/CMutationTCLib.hpp
@@ -58,6 +58,8 @@ class CMutationTCLib : public CNEMOGas {
 
   su2double Tref;                         /*!< \brief Reference temperature. */
 
+  bool NEWTON;                             /*!< \brief Boolean to use Newton-Raphson h/cp or Perturbation Method to calculate Jacobian. */
+
 public:
 
   /*!
diff --git a/SU2_CFD/src/fluid/CMutationTCLib.cpp b/SU2_CFD/src/fluid/CMutationTCLib.cpp
index 518ab7ce7e43..675427161816 100644
--- a/SU2_CFD/src/fluid/CMutationTCLib.cpp
+++ b/SU2_CFD/src/fluid/CMutationTCLib.cpp
@@ -56,6 +56,7 @@ CMutationTCLib::CMutationTCLib(const CConfig* config, unsigned short val_nDim):
   else if (Kind_TransCoeffModel == TRANSCOEFFMODEL::CHAPMANN_ENSKOG)
     transport_model = "Chapmann-Enskog_LDLT";
 
+  NEWTON = config->Get2TNewton();
   if (NoneqStateModel == "2T") {
       opt.setStateModel("ChemNonEqTTv");
   }
@@ -137,7 +138,7 @@ void CMutationTCLib::SetTDStateRhosTTv(vector<su2double>& val_rhos, su2double va
 
   Pressure = ComputePressure();
 
-  mix->setState(rhos.data(), temperatures.data(), 1);
+  mix->setState(rhos.data(), temperatures.data(), 1, NEWTON);
 
 }
 
@@ -295,7 +296,7 @@ vector<su2double>& CMutationTCLib::ComputeTemperatures(vector<su2double>& val_rh
   energies[0] = rhoE - rhoEvel;
   energies[1] = rhoEve;
 
-  mix->setState(rhos.data(), energies.data(), 0);
+  mix->setState(rhos.data(), energies.data(), 0, NEWTON);
 
   mix->getTemperatures(temperatures.data());
 
diff --git a/config_template.cfg b/config_template.cfg
index e00c16d41fbc..f556b44a694e 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -425,6 +425,8 @@ GAS_COMPOSITION= (0.77, 0.23, 0.0, 0.0, 0.0)
 FROZEN_MIXTURE= NO
 % 1-Temperature vs. 2-Temperature Model for Mutation++
 NONEQ_STATE_MODEL= 2T
+% Use Mutationpp h/cp relation for solving for 2T temperatures vs. Perturbation Method
+2T_NEWTON= YES
 %
 % Datadriven fluid model
 % For data-driven fluid models, an interpolation algorithm is used to retrieve the thermodynamic state for a given density and internal energy
diff --git a/subprojects/Mutationpp b/subprojects/Mutationpp
index c3650938b604..538bfbbcd032 160000
--- a/subprojects/Mutationpp
+++ b/subprojects/Mutationpp
@@ -1 +1 @@
-Subproject commit c3650938b604b0e1cdf84a7fe257eef10c7c5e29
+Subproject commit 538bfbbcd032b2059f8ad88ba2164c42db062463