diff --git a/src/problem/Formulation.cpp b/src/problem/Formulation.cpp
index 4a9cbe0263a56586502477330760acfeb7010dd6..c3742f01ca7c8e20c8cdb883eeffa6f177d55d75 100644
--- a/src/problem/Formulation.cpp
+++ b/src/problem/Formulation.cpp
@@ -70,6 +70,20 @@ namespace gmshddm
       }
     }
 
+    template< class T_Scalar >
+    typename Formulation< T_Scalar >::sourceStatus Formulation< T_Scalar >::indexStatusInShots(unsigned id) const {
+    for (unsigned i = 0; i < _activeShots.size(); ++i) {
+        if (_numberedPhysicalSourceTerms[i].count(id) > 0) {
+            if (_activeShots[i]) {
+                return sourceStatus::ACTIVE;
+            } else {
+                return sourceStatus::INACTIVE;
+            }
+        }
+    }
+    return sourceStatus::NOT_IN_NUMBERED_TERMS;
+}
+
     template< class T_Scalar >
     gmshfem::algebra::Vector< T_Scalar > Formulation< T_Scalar >::_sortDofValues(const gmshfem::algebra::Vector< T_Scalar > &vector, const unsigned int fieldTag)
     {
@@ -241,6 +255,12 @@ namespace gmshddm
       _physicalSourceTerms.insert(termTag);
     }
 
+    template< class T_Scalar >
+    void Formulation< T_Scalar >::numberedPhysicalSourceTerm(const unsigned int id, const unsigned int termTag)
+    {
+      _numberedPhysicalSourceTerms.at(id).insert(termTag);
+    }
+
     template< class T_Scalar >
     const gmshfem::function::ScalarFunction< T_Scalar > Formulation< T_Scalar >::artificialSource(const gmshfem::function::ScalarFunction< T_Scalar > &f)
     {
@@ -265,15 +285,18 @@ namespace gmshddm
       for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
         if(mpi::isItMySubdomain(idom)) {
           for(auto it = _volume[idom]->begin(); it != _volume[idom]->end(); ++it) {
-            if(_physicalSourceTerms.find((*it)->tag()) != _physicalSourceTerms.end()) {
-              if(_physicalCommutator) {
+            unsigned termTag = (*it)->tag();
+            auto termStatus = indexStatusInShots(termTag);
+            if(_physicalSourceTerms.count(termTag) > 0 || termStatus != sourceStatus::NOT_IN_NUMBERED_TERMS) {
+              if(_physicalCommutator && termStatus != sourceStatus::INACTIVE) {
                 (*it)->activate();
               }
               else {
+                // Disable if all physical terms are disabled OR if it's a disabled shot
                 (*it)->deactivate();
               }
             }
-            else if(_artificialSourceTerms.find((*it)->tag()) != _artificialSourceTerms.end()) {
+            else if(_artificialSourceTerms.count(termTag)) {
               if(_artificialCommutator) {
                 (*it)->activate();
               }
@@ -306,6 +329,37 @@ namespace gmshddm
       }
     }
 
+    template< typename T_Scalar >
+    void Formulation< T_Scalar >::setShotNumber(unsigned N)
+    {
+      _activeShots.resize(N, false);
+      _numberedPhysicalSourceTerms.resize(N);
+    }
+
+    template< typename T_Scalar >
+    void Formulation< T_Scalar >::enableShot(unsigned i)
+    {
+      _activeShots.at(i) = true;
+    }
+
+    template< typename T_Scalar >
+    void Formulation< T_Scalar >::disableShot(unsigned i)
+    {
+      _activeShots.at(i) = false;
+    }
+
+    template< typename T_Scalar >
+    void Formulation< T_Scalar >::toggleShot(unsigned i)
+    {
+      _activeShots.at(i) = !_activeShots.at(i);
+    }
+
+    template< typename T_Scalar >
+    bool Formulation< T_Scalar >::isShotEnabled(unsigned i) const
+    {
+      return _activeShots.at(i);
+    }
+
     template< class T_Scalar >
     void Formulation< T_Scalar >::addInterfaceField(field::InterfaceFieldInterface< T_Scalar > &interfaceField)
     {
diff --git a/src/problem/Formulation.h b/src/problem/Formulation.h
index ec02567d8ffdc6d29fc62e4cba5f01404ba20f27..23466c2f75eaa17bacc377ff225a98bb9eced66e 100644
--- a/src/problem/Formulation.h
+++ b/src/problem/Formulation.h
@@ -53,6 +53,18 @@ namespace gmshddm
       bool _IA;
       std::set< unsigned int > _physicalSourceTerms;
       std::set< unsigned int > _artificialSourceTerms;
+      std::vector< std::set< unsigned int > > _numberedPhysicalSourceTerms;
+      std::vector< bool > _activeShots;
+
+
+      enum class sourceStatus {
+        ACTIVE,
+        INACTIVE,
+        NOT_IN_NUMBERED_TERMS
+      };
+      sourceStatus indexStatusInShots(unsigned id) const;
+
+      
 
      protected:
       void _runIterationOperations(const unsigned int iteration, const gmshfem::scalar::Precision< T_Scalar > rnorm);
@@ -83,6 +95,15 @@ namespace gmshddm
       void artificialSourceTerm(const unsigned int termTag);
       void togglePhysicalAndArtificialSourceTerms();
 
+      // Multi source options
+      void setShotNumber(unsigned N);
+      void enableShot(unsigned i);
+      void disableShot(unsigned i);
+      void toggleShot(unsigned i);
+      bool isShotEnabled(unsigned i) const;
+      void numberedPhysicalSourceTerm(const unsigned int id, const unsigned int termTag);
+
+
       void addInterfaceField(field::InterfaceFieldInterface< T_Scalar > &interfaceField);
       void addInterfaceField(field::InterfaceFieldInterface< T_Scalar > &interfaceField, field::InterfaceFieldInterface< T_Scalar > &interfaceFieldMapping);
       field::InterfaceFieldInterface< T_Scalar > *getInterfaceField(const std::string &name) const;