diff --git a/docs/experimental-features.md b/docs/experimental-features.md
new file mode 100644
index 0000000..5a45287
--- /dev/null
+++ b/docs/experimental-features.md
@@ -0,0 +1,69 @@
+# Experimental Features
+
+## About the init system
+
+ElecSolver provides an init system feature that allows for users to obtain the inital condition of the electric system right when the simulation starts.
+
+This feature is essential when assessing the capacity of the system's reponse to a threshold input signal. But some configurations make the circuit not able to answer a threshold input:
+
+- A voltage source with capacitance connected to its terminals.
+- Inductances that, if disconnected from the circuit, split the current source terminals from one another
+
+In these cases, the init system is not defined and the system can only answer to differentiable inputs.
+
+When these conditions are not met the system can answer to a threshold input.
+
+When it is working the init system provides the solution at time $t=0+$. However, this feature does not always work:
+- It fails when having two capacitances in parallel in the system
+- It fails when having an inductance that is connected to a dead end in the circuit.
+- It fails when there are multiple inductances connected in series
+
+To circumvent these limitations we provided a workaround by providing a robust init system that can handle these special cases more effectively. Should your init system fail under the standard conditions, the robust init system can be used as an alternative.
+
+## About the robust init system
+
+The robust init system does not provide a physical solution at $t=0+$ in the same sense as the standard init system. Instead, it offers a mathematically consistent initial condition that allows the simulation to run even in cases where the standard init system would fail. It relies on the fact that currents in capacitors and potentials in terminals of lonely inductances are algebraic unknowns.
+
+It is important to note that while the robust init system increases the range of circuits that can be simulated, a user should not rely on its solution as a substitute for the physically accurate initial conditions provided by the standard init system.
+
+In a few words **do not take the result of the robust init system in your simulation as the first time step solution**. Try to use the standard init system whenever possible which is more physically sound.
+
+We are still trying to unify both solutions into a single coherent framework that can handle a wider range of initial conditions while maintaining physical accuracy.
+
+
+## Example Usage
+```python
+## Simple tetrahedron
+res_coords  = np.array([[1],[2]],dtype=int)
+res_data = np.array([1],dtype=float)
+## previous init system could not handle coils connecting a single node
+coil_coords  = np.array([[],[]],dtype=int)
+coil_data = np.array([],dtype=float)
+## previous init system could not handle 2 capacitors connected to the same node
+capa_coords = np.array([[0],[1]],dtype=int)
+capa_data = np.array([1],dtype=float)
+
+## total impedance
+mutual_coords=np.array([[],[]],dtype=int)
+mutual_data = np.array([],dtype=float)
+
+
+res_mutual_coords=np.array([[],[]],dtype=int)
+res_mutual_data = np.array([],dtype=float)
+
+elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutual_coords,mutual_data,res_mutual_coords,res_mutual_data)
+elec_sys.add_voltage_source(10,0,2)
+elec_sys.set_ground(0)
+elec_sys.build_system()
+S_i_robust,b = elec_sys.get_robust_init_system()
+sol = spsolve(S_i_robust.tocsr(),b)
+S_i, rhs = elec_sys.get_init_system()
+sol_init = spsolve(S_i.tocsr(),rhs)
+print("The two init system solutions are different")
+print("Robust init system solution:",sol)
+print("Standard init system solution:",sol_init)
+S1,S2,b = elec_sys.get_system()
+dt = 0.1
+print("Solving first timestep of the system with both init system solutions provide the same solution")
+np.testing.assert_allclose(spsolve(S2+dt*S1,b*dt+S2@sol), spsolve(S2+dt*S1,b*dt+S2@sol_init))
+```
diff --git a/docs/index.md b/docs/index.md
index 684cecb..e45c0c7 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -76,3 +76,4 @@ Then presents extra features:
 - [Solver suggestions](solver-suggestions.md)
 - [Extra uses: Hydraulic or Thermal system modeling](extra-uses-hydraulic-or-thermal-system-modeling.md)
 - [Netlist import feature](netlist-import-feature.md)
+- [Experimental features](experimental-features.md)
diff --git a/src/ElecSolver/FrequencySystemBuilder.py b/src/ElecSolver/FrequencySystemBuilder.py
index b3ce98d..203cb67 100644
--- a/src/ElecSolver/FrequencySystemBuilder.py
+++ b/src/ElecSolver/FrequencySystemBuilder.py
@@ -53,7 +53,7 @@ def graph_analysis(self):
         # Number of currents.
         self.number_intensities = self.all_impedences.shape[0]
         # Keep the connected subgraphs.
-        self.list_of_subgraphs = [ sub.tolist() for sub in compute_graph_components(self.all_coords)]
+        self.list_of_subgraphs = [ sub.tolist() for sub in compute_graph_components(self.all_coords,self.size)]
         self.number_of_subsystems = len(self.list_of_subgraphs)
         # Reassign previous grounds when repeating the analysis.
         if self.analysed:
diff --git a/src/ElecSolver/TemporalSystemBuilder.py b/src/ElecSolver/TemporalSystemBuilder.py
index 0353517..d969305 100644
--- a/src/ElecSolver/TemporalSystemBuilder.py
+++ b/src/ElecSolver/TemporalSystemBuilder.py
@@ -1,10 +1,10 @@
 import numpy as np
 from scipy.sparse import coo_matrix, block_diag, coo_array
-from .utils import SolutionTemporal,GradientsParametersTemporal,compute_graph_components
+from .utils import SolutionTemporal,GradientsParametersTemporal,compute_graph_components,get_roots
 import warnings
 
 class TemporalSystemBuilder():
-    def __init__(self,coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,inductive_mutual_coords,inductive_mutual_data,res_mutual_coords,res_mutual_data):
+    def __init__(self,coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,inductive_mutual_coords,inductive_mutual_data,res_mutual_coords,res_mutual_data,size=None):
         """Initialize a temporal electrical system builder.
 
         Parameters
@@ -35,6 +35,8 @@ def __init__(self,coil_coords,coil_data,res_coords,res_data,capa_coords,capa_dat
         res_mutual_data : numpy.ndarray, shape (P,)
             Mutual resistance values. Each sign follows the node order in
             ``coil_coords`` or ``res_coords``.
+        size : int, optional
+            The size of the system. If not provided, it will be determined during graph analysis.
         """
 
         self.coil_coords=coil_coords
@@ -58,6 +60,7 @@ def __init__(self,coil_coords,coil_data,res_coords,res_data,capa_coords,capa_dat
         self.rhs = (np.array([]),(np.array([],dtype=int),))
 
         self.analysed=False
+        self.size = size
 
     def graph_analysis(self):
         """Analyze connected components and initialize system dimensions.
@@ -84,11 +87,14 @@ def graph_analysis(self):
         self.all_impedences = np.concatenate([self.coil_data,self.res_data,self.capa_data,self.voltage_source_data],axis=0)
 
         # Actual number of nodes in the system.
-        self.size = np.max(self.all_coords)+1
+        if self.size is None:
+            self.size = np.max(self.all_coords)+1
+        else:
+            self.size = max(np.max(self.all_coords)+1, self.size)
         # Number of currents.
         self.number_intensities = self.all_impedences.shape[0]
         # Keep the connected subgraphs.
-        self.list_of_subgraphs = [ list(sub) for sub in compute_graph_components(self.all_coords)]
+        self.list_of_subgraphs = [ list(sub) for sub in compute_graph_components(self.all_coords,self.size)]
         self.number_of_subsystems = len(self.list_of_subgraphs)
         # Reassign previous grounds when repeating the analysis.
         if self.analysed:
@@ -455,6 +461,146 @@ def get_init_system(self,sparse_rhs=True):
             np.add.at(rhs, nodes, data_rhs)
         return sys,rhs
 
+    def get_robust_init_system(self,sparse_rhs=True):
+        """Fonction to return a more robust initial system and the right hand side
+        in case the initial system is singular, this comes at the cost of additionnal
+        computational cost. Compared to the standart init system the initial condition it provided
+        is at time t=0^{-} whereas the initial system assumes t=0^{+}.
+        Therefore when using the robust initial system one should be aware that it does not provide a sound initial condition at t=0^{+}.
+
+        Basically to make the it work it needs:
+            - to merge the capacitors nodes
+            - to remove the coils wires, the mutuals for the coils
+            - to make sure the merged nodes have same potential (really needed?)
+        Returns
+        -------
+        sys : scipy.sparse.coo_matrix, shape (n, n)
+            Robust initial-condition system matrix.
+        rhs : scipy.sparse.coo_array or numpy.ndarray, shape (n,)
+            Right-hand side of the robust initial system.
+        """
+        # Run graph analysis if needed.
+        if not self.analysed:
+            self.graph_analysis()
+        # Assign grounds if needed.
+        self.affect_potentials()
+        ## First compute roots for the capacitance graph
+        capa_roots = get_roots(self.capa_coords,self.size)
+        reduced_tension_coords = capa_roots[self.voltage_source_coords]
+        if any(reduced_tension_coords[0]==reduced_tension_coords[1]):
+            for i,(a,b) in enumerate(zip(reduced_tension_coords[0],reduced_tension_coords[1])):
+                if a==b:
+                    print(f"Voltage source {i} is bridged by capacitors.")
+            raise ValueError("System is not valid for threshold input.")
+
+
+        ## Second computing graph for res, capa + tension_sources
+        neg_induct_roots = get_roots(np.concatenate([self.res_coords,self.capa_coords,self.voltage_source_coords],axis=1),self.size)
+        reduced_current_coords = neg_induct_roots[self.current_source_coords]
+        if any(reduced_current_coords[0]!=reduced_current_coords[1]):
+            for i,(a,b) in enumerate(zip(reduced_current_coords[0],reduced_current_coords[1])):
+                if a!=b:
+                    print(f"Current source {i}'s input and output don't belong to the same subsystem or are split by an inductance.")
+            raise ValueError("System is not valid for threshold input.")
+
+        ## Merging nodes that belong to the same capacitance graph
+        ## Resistances
+        local_res_coords = capa_roots[self.res_coords]
+        res_to_delete = local_res_coords[0]==local_res_coords[1]
+        local_res_coords = local_res_coords[:,~res_to_delete]
+        local_res_data = self.res_data[~res_to_delete]
+        # Need to keep track of deleted resistances for building back the augmented system later
+        rescaler_res = np.cumsum(res_to_delete.astype(int))[~res_to_delete]
+
+        ## Current sources
+        local_current_source_coords = capa_roots[self.current_source_coords]
+        current_source_to_delete = local_current_source_coords[0]==local_current_source_coords[1]
+        local_current_source_coords = self.current_source_coords[:,~current_source_to_delete]
+        local_current_source_data = self.current_source_data[~current_source_to_delete]
+
+        ## Voltage sources
+        local_voltage_source_coords = capa_roots[self.voltage_source_coords]
+        local_voltage_source_data = self.voltage_source_data
+
+
+        builder = TemporalSystemBuilder(
+            coil_coords=np.array([[],[]],dtype=int),
+            coil_data=np.array([],dtype=float),
+            res_coords=local_res_coords,
+            res_data=local_res_data,
+            capa_coords=np.array([[],[]],dtype=int),
+            capa_data=np.array([],dtype=float),
+            inductive_mutual_coords=np.array([[],[]],dtype=int),
+            inductive_mutual_data=np.array([],dtype=float),
+            res_mutual_coords=np.array([[],[]],dtype=int),
+            res_mutual_data=np.array([],dtype=float),
+            size=self.size
+        )
+        ## setting the sources manually
+        builder.current_source_coords = local_current_source_coords
+        builder.current_source_data = local_current_source_data
+        builder.voltage_source_coords = local_voltage_source_coords
+        builder.voltage_source_data = local_voltage_source_data
+        builder.source_count=local_voltage_source_data.shape[0]
+        builder.set_ground(*capa_roots[self.affected_potentials])
+        ## building the condensed system: merged capacitance, removed coils, removed short-cut resistances,
+        ## removed short cut current sources
+        builder.build_system()
+        S_init_condensed = builder.S1
+        rhs = builder.rhs
+        ## expanding the condensed system to the full system
+        rescaler_res = np.arange(self.coil_data.shape[0],self.coil_data.shape[0]+local_res_data.shape[0]) +rescaler_res
+        rescaler_voltage_source = np.arange(self.coil_data.shape[0]+self.res_data.shape[0]+ self.capa_data.shape[0],self.coil_data.shape[0]+self.res_data.shape[0]+ self.capa_data.shape[0]+self.voltage_source_data.shape[0])
+        rescaler_potentials = self.coil_data.shape[0]+self.res_data.shape[0]+self.capa_data.shape[0]+self.voltage_source_data.shape[0] + np.arange(0,self.size)
+        j_rescaler = np.concatenate([rescaler_res, rescaler_voltage_source, rescaler_potentials])
+        data,(i_s,j_s) = S_init_condensed
+        data_rhs,(i_rhs,) = rhs
+        j_s = j_rescaler[j_s]
+        equation_index = i_s.max()+1
+        ## cancel currents in coils, cancel currents in capacities and cancel in short cut resistances
+        j_s_coils = np.arange(self.coil_data.shape[0])
+        i_s_coils = equation_index + np.arange(self.coil_data.shape[0])
+        data_coils = np.ones(self.coil_data.shape[0],dtype=float)
+        ## cancel currents in short-cut
+        j_s_res = (np.arange(self.res_data.shape[0])+self.coil_data.shape[0])[res_to_delete]
+        i_s_res = equation_index + i_s_coils.shape[0] + np.arange(j_s_res.shape[0])
+        data_res = np.ones(j_s_res.shape[0],dtype=float)
+        ## cancel current in capacitances
+        j_s_capa = (np.arange(self.capa_data.shape[0])+self.coil_data.shape[0]+self.res_data.shape[0])
+        i_s_capa = equation_index +i_s_coils.shape[0]+ i_s_res.shape[0] + np.arange(j_s_capa.shape[0])
+        data_capa = np.ones(j_s_capa.shape[0],dtype=float)
+
+        affected_potentials_reduced = np.array(builder.affected_potentials,dtype=int)
+        mask = capa_roots[affected_potentials_reduced] != affected_potentials_reduced
+        i_s_potential = builder.size + builder.number_intensities - builder.source_count -builder.number_of_subsystems + np.arange(affected_potentials_reduced.shape[0])
+        j_s_potential = j_rescaler[builder.offset_j + capa_roots[affected_potentials_reduced]]
+        data_potential = -np.ones(affected_potentials_reduced.shape[0],dtype=float)
+
+        i_s_potential = i_s_potential[mask]
+        j_s_potential = j_s_potential[mask]
+        data_potential = data_potential[mask]
+
+        ## Add all equations
+        j_s = np.concatenate([j_s,j_s_coils, j_s_res, j_s_capa, j_s_potential])
+        i_s = np.concatenate([i_s,i_s_coils, i_s_res, i_s_capa, i_s_potential])
+        data = np.concatenate([data,data_coils, data_res, data_capa, data_potential])
+
+        j_s = j_s.astype(int)
+        i_s = i_s.astype(int)
+
+        size = self.number_intensities+self.size
+        sys = coo_matrix((data,(i_s,j_s)),shape=(size,size))
+        if sparse_rhs:
+            rhs = coo_array((data_rhs,(i_rhs,)),shape=(size,))
+            rhs.sum_duplicates()
+        else:
+            rhs = np.zeros(size)
+            (data_rhs,(i_rhs,)) = self.rhs
+            np.add.at(rhs, i_rhs, data_rhs)
+        ## rhs here is only for the init system do not mix up with the main system's rhs
+        return sys,rhs
+
+
     def get_system(self,sparse_rhs=True):
         """Return the real and derivative systems and right-hand side.
 
diff --git a/src/ElecSolver/utils.py b/src/ElecSolver/utils.py
index e5331a0..7354d7e 100644
--- a/src/ElecSolver/utils.py
+++ b/src/ElecSolver/utils.py
@@ -1,5 +1,6 @@
 import numpy as np
 from scipy.sparse import coo_matrix
+from scipy.sparse.csgraph import connected_components
 from collections import namedtuple
 
 SolutionFrequency = namedtuple("ComplexSolution",["intensities","potentials","intensities_sources"])
@@ -127,37 +128,50 @@ def constant_block_diag(A,repetitions):
     data = np.tile(A.data,(repetitions,))
     return coo_matrix((data,(rows,cols)),shape=(repetitions*size,repetitions*size))
 
-def compute_graph_components(all_coords):
+def get_roots(all_coords, number_of_nodes=None):
+    """Find the root of each node in a graph.
+
+    Parameters
+    ----------
+    all_coords : numpy.ndarray, shape (2, n_edges)
+        Node indices for each graph edge.
+    number_of_nodes : int, optional
+        Total number of nodes in the graph. Defaults to ``max(all_coords) + 1``.
+
+    Returns
+    -------
+    roots : numpy.ndarray, shape (number_of_nodes,)
+        Root of each node: the lowest node index of its connected component.
+    """
+    if number_of_nodes is None:
+        number_of_nodes = int(np.max(all_coords))+1
+    number_of_edges = all_coords.shape[1]
+    graph = coo_matrix((np.ones(number_of_edges,dtype=np.int8),(all_coords[0],all_coords[1])),
+                       shape=(number_of_nodes,number_of_nodes))
+    _, labels = connected_components(graph,directed=False)
+    # Labels are assigned in node order, so the first occurrence is the lowest index.
+    _, first = np.unique(labels,return_index=True)
+    return first[labels]
+
+def compute_graph_components(all_coords,number_of_nodes=None):
     """Compute graph components from edge coordinates.
+    It is a simple wrapper around :func:`get_roots` to group nodes by their connected components.
 
     Parameters
     ----------
     all_coords : numpy.ndarray, shape (2, n_edges)
         Node indices for each graph edge.
+    number_of_nodes : int, optional
+        Total number of nodes in the graph. Defaults to ``max(all_coords) + 1``.
 
     Returns
     -------
     list of numpy.ndarray
         Node indices in each connected component.
     """
-    number_of_nodes = np.max(all_coords)+1
-    number_of_edges = all_coords.shape[1]
-    # Initialize the root of each node to itself.
-    roots = np.arange(0,number_of_nodes,dtype=int)
-    # Make each point refer to its lowest-index neighbor.
-    for i in range(number_of_edges):
-        min_value = min(all_coords[0,i],all_coords[1,i])
-        roots[all_coords[0,i]] = min(min_value,roots[all_coords[0,i]])
-        roots[all_coords[1,i]] = min(min_value,roots[all_coords[1,i]])
-
-    converged = False
-    # Iteratively update each node to point to the root of its component.
-    while not converged:
-        previous_roots = roots.copy()
-        roots = roots[roots]
-        # Converged when no roots have changed in this iteration.
-        converged = np.all(previous_roots==roots)
-
+    if number_of_nodes is None:
+        number_of_nodes = int(np.max(all_coords))+1
+    roots = get_roots(all_coords,number_of_nodes)
     return [np.arange(0,number_of_nodes,dtype=int)[roots==i] for i in np.unique(roots)]
 
 
diff --git a/tests/test_temporal_system.py b/tests/test_temporal_system.py
index a0277ef..597e801 100644
--- a/tests/test_temporal_system.py
+++ b/tests/test_temporal_system.py
@@ -148,7 +148,7 @@ def test_one_shot_temporal():
     elec_sys.set_ground(0)
     elec_sys.build_system()
 
-    S_i,b = elec_sys.get_init_system(sparse_rhs=False)
+    S_i,b = elec_sys.get_robust_init_system(sparse_rhs=False)
     sol = spsolve(S_i.tocsr(),b)
     S1,S2,rhs = elec_sys.get_system(sparse_rhs=False)
     dt=0.08
@@ -436,6 +436,64 @@ def test_reverting_named_tuple():
         named_sol = elec_sys.build_intensity_and_voltage_from_vector(sol)
         sol = elec_sys.build_vector_from_intensity_and_voltage(named_sol)
 
+def test_robust_init_system():
+    ## Simple tetrahedron
+    res_coords  = np.array([[2,0],[3,1]],dtype=int)
+    res_data = np.array([1,1],dtype=float)
+    ## previous init system could not handle coils connecting a single node
+    coil_coords  = np.array([[1,3],[2,5]],dtype=int)
+    coil_data = np.array([1,4],dtype=float)
+    ## previous init system could not handle 2 capacitors connected to the same node
+    capa_coords = np.array([[0,0],[1,1]],dtype=int)
+    capa_data = np.array([1,3],dtype=float)
+
+    ## total impedance
+    mutual_coords=np.array([[],[]],dtype=int)
+    mutual_data = np.array([],dtype=float)
+
+
+    res_mutual_coords=np.array([[],[]],dtype=int)
+    res_mutual_data = np.array([],dtype=float)
+
+    elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutual_coords,mutual_data,res_mutual_coords,res_mutual_data)
+    elec_sys.add_voltage_source(10,0,3)
+    elec_sys.set_ground(0,4)
+    elec_sys.build_system()
+    S_i,b = elec_sys.get_robust_init_system()
+    sol = spsolve(S_i.tocsr(),b)
+
+
+def test_compare_robust_init_system():
+    ## Simple tetrahedron
+    res_coords  = np.array([[1],[2]],dtype=int)
+    res_data = np.array([1],dtype=float)
+    ## previous init system could not handle coils connecting a single node
+    coil_coords  = np.array([[],[]],dtype=int)
+    coil_data = np.array([],dtype=float)
+    ## previous init system could not handle 2 capacitors connected to the same node
+    capa_coords = np.array([[0],[1]],dtype=int)
+    capa_data = np.array([1],dtype=float)
+
+    ## total impedance
+    mutual_coords=np.array([[],[]],dtype=int)
+    mutual_data = np.array([],dtype=float)
+
+
+    res_mutual_coords=np.array([[],[]],dtype=int)
+    res_mutual_data = np.array([],dtype=float)
+
+    elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutual_coords,mutual_data,res_mutual_coords,res_mutual_data)
+    elec_sys.add_voltage_source(10,0,2)
+    elec_sys.set_ground(0)
+    elec_sys.build_system()
+    S_i_robust,b = elec_sys.get_robust_init_system()
+    sol = spsolve(S_i_robust.tocsr(),b)
+    S_i, rhs = elec_sys.get_init_system()
+    sol_init = spsolve(S_i.tocsr(),rhs)
+    S1,S2,b = elec_sys.get_system()
+    dt = 0.1
+    np.testing.assert_allclose(spsolve(S2+dt*S1,b*dt+S2@sol), spsolve(S2+dt*S1,b*dt+S2@sol_init))
+
 
 if __name__ == "__main__":
     test_lonely_nodes()
@@ -446,3 +504,6 @@ def test_reverting_named_tuple():
     test_voltage()
     test_temporal()
     test_temporal2()
+    test_reverting_named_tuple()
+    test_robust_init_system()
+    compare_robust_init_system()
diff --git a/tests/test_utils.py b/tests/test_utils.py
index dd11bcc..ef35ec1 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -7,8 +7,13 @@ def test_block_diag():
     constant_block_diag(A,3).todense()
 
 def test_graph_components():
-    all_coords = np.array([[0,1],[1,2],[3,4]]).T
+    all_coords = np.array([[0,1,3],[1,2,4]],dtype=int)
     components = compute_graph_components(all_coords)
+    print(components)
     assert len(components)==2
     assert np.all(components[0]==np.array([0,1,2]))
-    assert np.all(components[1]==np.array([3,4]))
\ No newline at end of file
+    assert np.all(components[1]==np.array([3,4]))
+
+if __name__ == "__main__":
+    test_block_diag()
+    test_graph_components()
\ No newline at end of file
diff --git a/zensical.toml b/zensical.toml
index 6ccc4d0..4e80ea6 100644
--- a/zensical.toml
+++ b/zensical.toml
@@ -16,6 +16,7 @@ nav = [
   "solver-suggestions.md",
   "extra-uses-hydraulic-or-thermal-system-modeling.md",
   "netlist-import-feature.md",
+  "experimental-features.md",
 ]
 extra_javascript = [
     "javascripts/mathjax.js",