diff --git a/README.md b/README.md
index 2028084..3c7e004 100644
--- a/README.md
+++ b/README.md
@@ -10,16 +10,24 @@ This repository is **not** a general-purpose electrical system solver. Instead,
 - The graph-based description of an electric network
 - The corresponding sparse linear system to solve
 
+In a very simple way, ElecSolver takes as an input the Resistances $R$, Capacitances $C$, Inductances $L$, Mutuals $M$, Current sources $I$  and Voltage sources $V$ along with the connectivity graph $G$ of the system and outputs the linear system: matrix $S$ and vector $b$ to solve in order to get the solution of the electric problem.
+
+$$f_{\text{ElecSolver}}(R,C,L,M,I,V,G) = (S,b)$$
+
+It is then up to the user to choose the solver they want for solving the system $Sx=b$ and to manage the time iterations if needed for temporal problems.
+
+The main goal of ElecSolver is to provide a friendly Python interface for simulating and optimizing analog electric systems. While suitable for small circuit simulations, its strength lies in its scalability: it is able to build linear systems with millions of nodes and components.
+
 
-Its main goal is to provide a friendly Python interface for simulating analog electric systems. While suitable for small circuit simulations, its strength lies in its scalability: it is able to build linear systems with millions of nodes and components.
 
 
 > [!IMPORTANT]
 > ElecSolver has been designed with the following specifications in mind:
 > - The time needed for building the linear system must be negligible compared to the time needed for solving it.
-> - Handle natively inductive mutuals and resistive mutuals
-> - Handle as many coupled electric systems that one wants.
-> - Deal with lonely nodes and lonely edges in the electric graph: the problem can be well posed and thus solved.
+> - Handle natively inductive mutuals and resistive mutuals (skin effects).
+> - Handle as many coupled electric systems as one wants.
+> - Deal with lonely nodes and lonely edges in the electric graph: the problem can be well-posed and thus solved.
+> - Allow backpropagation of gradients through the system for optimization purposes.
 
 
 
@@ -27,12 +35,15 @@ Its main goal is to provide a friendly Python interface for simulating analog el
 > Non-linear components are not supported. You must manage event detection and system updates yourself.
 
 
-## Table of content
+## Table of contents
 
 - [ElecSolver](#elecsolver)
   - [Overview](#overview)
-  - [Table of content](#table-of-content)
+  - [Table of contents](#table-of-contents)
   - [How to install](#how-to-install)
+    - [Using pip](#using-pip)
+    - [Using conda/mamba](#using-condamamba)
+    - [From source](#from-source)
   - [Components](#components)
     - [FrequencySystemBuilder](#frequencysystembuilder)
       - [Features](#features)
@@ -44,20 +55,33 @@ Its main goal is to provide a friendly Python interface for simulating analog el
   - [Solver suggestions](#solver-suggestions)
   - [Extra uses: Hydraulic or Thermal system modeling](#extra-uses-hydraulic-or-thermal-system-modeling)
   - [Netlist import feature](#netlist-import-feature)
+  - [Gradient backpropagation](#gradient-backpropagation)
+    - [Example of backpropagation of TemporalSystemBuilder](#example-of-backpropagation-of-temporalsystembuilder)
+    - [Example of backpropagation of FrequencySystemBuilder](#example-of-backpropagation-of-frequencysystembuilder)
+    - [Conclusion on gradient backpropagation](#conclusion-on-gradient-backpropagation)
+
 
 ## How to install
-For now this package is distributed on pypi and can be installed using pip and conda/mamba
+### Using pip
+ElecSolver is distributed on PyPI and can be installed using pip or conda/mamba.
 ```
 pip install ElecSolver
 ```
-or
+### Using conda/mamba
+ElecSolver is also available on conda-forge and can be installed using conda or mamba:
 ```
 conda install elecsolver
 ```
-For solving the linear systems we advise using MUMPS through `python-mumps` available for linux, macOS and Windows. It can be installed via conda:
+For solving linear systems, we advise using MUMPS through `python-mumps`, available for Linux, macOS, and Windows. It can be installed via conda:
 ```
 conda install python-mumps
 ```
+### From source
+You can also install ElecSolver from source by cloning the repository and running:
+```
+pip install .
+```
+you will need to have `numpy`, `scipy` and `networkx` installed in your environment
 
 
 
@@ -74,6 +98,7 @@ This class handles **frequency-domain** analysis of linear electric systems.
 - Detects and couples multiple subsystems
 - Accepts arbitrary complex impedances and mutuals
 - Constructs sparse linear systems (COO format)
+- Allows backpropagation of gradients from the constructed linear system (`S` and `rhs`) to the parameters of the system for optimization purposes
 
 
 > [!TIP]
@@ -84,7 +109,7 @@ This class handles **frequency-domain** analysis of linear electric systems.
 We would like to study the following system:
 ![Multiple system](docs/img/schema.png)
 
-this can simply be defined in the following manner (We took R=1, L=1 and M=2):
+This can simply be defined in the following manner (we take $R=1$, $L=1$, and $M=2$):
 ```python
 import numpy as np
 from scipy.sparse.linalg import spsolve
@@ -112,7 +137,7 @@ electric_sys = FrequencySystemBuilder(
 # Add source (Current source here)
 electric_sys.add_current_source(intensity=10, input_node=2, output_node=0)
 # Set ground
-# 2 values because one for each subsystem
+# Two values because one is needed for each subsystem
 electric_sys.set_ground(0, 3)
 # Building system
 electric_sys.build_system()
@@ -127,10 +152,10 @@ frequencial_response = electric_sys.build_intensity_and_voltage_from_vector(sol)
 print(frequencial_response.potentials[3]-frequencial_response.potentials[4])
 ```
 #### Adding a Parallel Resistance
-We want to add components in parallel with existing components for instance inserting a resistor in parallel with the first inductance (between nodes 0 and 2)
+We want to add components in parallel with existing components, for instance by inserting a resistor in parallel with the first inductance (between nodes 0 and 2).
 ![Parallel system](docs/img/schema3.png)
 
-In python, simply add the resistance to the list of impedence in the very first lines of the script:
+In Python, simply add the resistance to the list of impedence in the very first lines of the script:
 
 ```python
 import numpy as np
@@ -138,7 +163,7 @@ from scipy.sparse.linalg import spsolve
 from ElecSolver import FrequencySystemBuilder
 
 
-# We add an additionnal resistance between 0 and 2
+# We add an additional resistance between 0 and 2
 impedence_coords = np.array([[0,0,1,3,0],[1,2,2,4,2]], dtype=int)
 impedence_data = np.array([1, 1j,1, 1j,1], dtype=complex)
 
@@ -160,6 +185,7 @@ This class models **time-dependent** systems using resistors, capacitors, coils,
 - Detects and couples multiple subsystems
 - Accepts 3 dipole types: resistances, capacities and coils
 - Constructs sparse linear systems (COO format)
+- Allows backpropagation of gradients from the constructed linear systems (plural! gradients from `S_init`, `S1`, `S2` and `rhs`) to the parameters of the system for optimization purposes
 
 #### Example
 
@@ -167,7 +193,7 @@ This class models **time-dependent** systems using resistors, capacitors, coils,
 We would like to study the following system:
 ![Temporal system](docs/img/schema2.png)
 
-with R=1, L=0.1, C=2 this gives:
+With $R=1$, $L=0.1$, and $C=2$, this gives:
 
 ```python
 import numpy as np
@@ -243,13 +269,13 @@ This outputs the following graph that displays the intensity passing through the
 
 ## Extra uses: Hydraulic or Thermal system modeling
 
-This repository can be used as is in order to model the mass flow or thermal flux in respectively Hydraulic networks or Thermal networks where a difference of pressure or a difference of temperature can be assimilated to a tension source. Since electric potentials are always computed relatively to the ground node you might need to rescale the resulting potentials:
+This repository can be used as is to model mass flow or thermal flux in hydraulic or thermal networks, respectively, where a pressure difference or a temperature difference can be assimilated to a tension source. Since electric potentials are always computed relative to the ground node, you may need to rescale the resulting potentials:
 
 We are considering the following hydraulic problem:
 
 ![Hydraulic system](docs/img/hydraulic.png)
 
-Taking R=1 this gives
+Taking $R=1$ this gives
 
 ```python
 import numpy as np
@@ -304,7 +330,7 @@ print("Debit through the system",solution.intensities_sources[0])
 ```
 ## Netlist import feature
 
-A new class, named NetlistParser allows importing passive netlist and building a TemporalSystem instance.
+A new class named NetlistParser allows importing a passive netlist and building a TemporalSystemBuilder instance.
 Solving the system can then be performed like any other example above.
 
 ```python
@@ -376,4 +402,130 @@ plt.plot(arange(steps, dtype=float)*dt, voltage_src, label="V(1-0)") # equal to
 plt.legend()
 plt.savefig("intensities_res.png")
 ```
+## Gradient backpropagation
+
+As we said in the introduction, ElecSolver is roughly a function that takes as an input the graph-based description of an electric system and outputs the linear system to solve:
+
+$$f_{\text{ElecSolver}}(R,C,L,M,I,V,G) = (S,b)$$
+
+
+In order to allow optimization of the system, ElecSolver now gives the possibility to backpropagate gradients from the linear system to the parameters of the system (i.e., by giving $\frac{\partial R}{\partial S}$ , $\frac{\partial R}{\partial b}$ , $\frac{\partial C}{\partial S}$ , $\frac{\partial C}{\partial b}$ , etc.). This allows one to use ElecSolver in an optimization loop where one wants to optimize the parameters of the system with gradient descent, for instance.
+
+### Example of backpropagation of TemporalSystemBuilder
+In this example, we want to optimize the capacity values of the system in order to reach a measured target current for `t=0.8`.
+We can do this by backpropagating the gradients from the solution of the linear system to the capa_data array of the TemporalSystemBuilder instance and then performing gradient descent on capa_data.
+```python
+import numpy as np
+from scipy.sparse.linalg import spsolve
+from ElecSolver import TemporalSystemBuilder
+## Simple tetrahedron
+res_coords  = np.array([[0,2],[1,3]],dtype=int)
+res_data = np.array([1,1],dtype=float)
+
+coil_coords  = np.array([[1,0],[2,3]],dtype=int)
+coil_data = np.array([1,1],dtype=float)
+
+capa_coords = np.array([[1,2],[3,0]],dtype=int)
+capa_data = np.array([1,1],dtype=float)
+## The target solution we want to reach is the solution of the system when capa_data = np.array([0.1,1],dtype=float)
+
+## mutuals
+mutuals_coords=np.array([[],[]],dtype=int)
+mutuals_data = np.array([],dtype=float)
+
+
+res_mutuals_coords=np.array([[],[]],dtype=int)
+res_mutuals_data = np.array([],dtype=float)
+
+## initializing system
+elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+elec_sys.add_current_source(10,1,0)
+elec_sys.set_ground(0)
+elec_sys.build_system()
+
+# Getting initial condition system
+S_i,rhs = elec_sys.get_init_system(sparse_rhs=True)
+# Getting temporal systems
+S1,S2,rhs = elec_sys.get_system(sparse_rhs=True)
+# Solving initial condition system
+sol_init =spsolve(S_i.tocsr(),rhs.todense())
+## Time iteration with euler implicit scheme for 1 timestep
+dt=0.8
+B = rhs*dt+S2@sol_init
+A = S2+dt*S1
+sol = spsolve(A,B)
+## Target solution (artificially made by setting capa_data = np.array([0.1,1],dtype=float))
+sol_target = [ 3.24786325, -1.1965812,  -6.16809117,  0.61253561,  0.58404558, -2.63532764, 0., 6.16809117,  2.10826211,  1.4957265 ]
+for i in range(1000):
+
+    ## computing gradients
+    dB = 2*spsolve(A.T, sol - sol_target)
+    ## chain rule for gradients of capa_data (S2 appears twice in the computation graph)
+    dS2 = -( dB[S2.row]*sol[S2.col]) ## Gradient from the solver (A\B) to S2
+    dS2 += (dB[S2.row]*sol_init[S2.col]) ## Gradient from the solver (A\B) to S2 through B
+
+    ## Backpropagate gradients from dS2 to capa_data
+    gradients = elec_sys.backpropagate_gradients(dS2=dS2)
+
+    ## change the values of capa_data using gradient descent
+    elec_sys.capa_data = elec_sys.capa_data - 0.01*gradients.capa_data
+    ## After the update of capa_data we need to rebuild the system to update S1, S2 and rhs
+    elec_sys.build_system()
+    S1,S2,rhs = elec_sys.get_system(sparse_rhs=True)
+    ## recomputing solution with euler implicit scheme
+    B = rhs*dt+S2@sol_init
+    A = S2+dt*S1
+    sol = spsolve(A,B)
+## Checking whether we converged to the right solution
+np.testing.assert_allclose(elec_sys.capa_data, np.array([0.1,1],dtype=float))
+```
+
+The function `backpropagate_gradients` of TemporalSystemBuilder allows to backpropagate gradients from `S_init`, `S1`, `S2` and `rhs`. See tests.test_gradients for more examples of backpropagation of gradients from different systems.
+
+### Example of backpropagation of FrequencySystemBuilder
+We want to optimize the value of a tension source in order to reach a target solution for the frequency response of the system.
+```python
+import numpy as np
+from scipy.sparse.linalg import spsolve
+from ElecSolver import FrequencySystemBuilder
+
+## sparse python res matrix
+impedence_coords = np.array([[0,0,1],[1,2,2]],dtype=int)
+impedence_data = np.array([1,1,1],dtype=complex)
+
+## mutuals
+mutuals_coords=np.array([[0],[1]],dtype=int)
+mutuals_data = np.array([2.j],dtype=complex)
+
+
+electric_sys = FrequencySystemBuilder(impedence_coords,impedence_data,mutuals_coords,mutuals_data)
+## target solution is the solution of the system when voltage=5
+electric_sys.add_voltage_source(voltage=10,input_node=1,output_node=0)
+# setting the ground
+electric_sys.set_ground(0)
+electric_sys.build_system()
+
+## getting system
+sys,b = electric_sys.get_system(sparse_rhs=True)
+sol = spsolve(sys.tocsr(),b.todense())
+## Target solution (Artificially made by setting voltage_source_data = np.array([5],dtype=complex))
+sol_target = np.array([-1.66666667+1.66666667j, -0.83333333+1.66666667j, 0.83333333-1.66666667j, 2.5-3.33333333j, 0.+0.j, 5+0.j, 4.16666667+1.66666667j])
+for i in range(3000):
+    ## computing gradients of the squared error between sol and sol_target with respect to the second member b of the system (the one that contains the voltage source data)
+    db = 2*spsolve(sys.tocsr().conj().T, sol - sol_target)
+    drhs = db[b.row]
+    ## Backpropagate gradients from drhs to voltage_source_data
+    gradients = electric_sys.backpropagate_gradients(drhs=drhs)
+    ## Performing gradient descent on voltage_source_data
+    electric_sys.voltage_source_data = electric_sys.voltage_source_data - 0.01*gradients.voltage_source_data
+    ## After the update of voltage_source_data we need to rebuild the system to update sys and b
+    electric_sys.build_system()
+    sys,b = electric_sys.get_system(sparse_rhs=True)
+    sol = spsolve(sys.tocsr(),b.todense())
+## Checking whether we converged to the right solution
+np.testing.assert_allclose(electric_sys.voltage_source_data, np.array([5],dtype=complex))
+```
+The function `backpropagate_gradients` of FrequencySystemBuilder allows to backpropagate gradients from `S`, and `rhs`. See tests.test_gradients for more examples of backpropagation of gradients from different systems.
 
+### Conclusion on gradient backpropagation
+Although it provides the backpropagation feature, ElecSolver does not provide an automatic differentiation mechanism. Many libraries can be used to this extent, such as `autograd`, `jax`, `PyTorch` and many more, to avoid the hassle of computing gradients manually. ElecSolver aims to limit its dependencies to keep its versatility and remain a simple, user-friendly tool.
diff --git a/docs/components/frequency-system-builder.md b/docs/components/frequency-system-builder.md
index b7f5572..0146318 100644
--- a/docs/components/frequency-system-builder.md
+++ b/docs/components/frequency-system-builder.md
@@ -87,3 +87,76 @@ impedence_data = np.array([1, 1j, 1, 1j, 1], dtype=complex)
 mutuals_coords = np.array([[1], [3]], dtype=int)
 mutuals_data = np.array([2.0j], dtype=complex)
 ```
+
+## Gradient Backpropagation
+
+`FrequencySystemBuilder` can backpropagate gradients from `S` and `rhs` to model parameters.
+
+This enables gradient-based optimization loops directly on source values or component data.
+
+### Example: Optimize a Voltage Source
+
+In this example, we optimize the voltage source value so the system response matches a target solution.
+
+```python
+import numpy as np
+from scipy.sparse.linalg import spsolve
+from ElecSolver import FrequencySystemBuilder
+
+## sparse python res matrix
+impedence_coords = np.array([[0, 0, 1], [1, 2, 2]], dtype=int)
+impedence_data = np.array([1, 1, 1], dtype=complex)
+
+## mutuals
+mutuals_coords = np.array([[0], [1]], dtype=int)
+mutuals_data = np.array([2.0j], dtype=complex)
+
+electric_sys = FrequencySystemBuilder(
+    impedence_coords,
+    impedence_data,
+    mutuals_coords,
+    mutuals_data,
+)
+
+# Target solution is the solution of the system when voltage=5
+electric_sys.add_voltage_source(voltage=10, input_node=1, output_node=0)
+electric_sys.set_ground(0)
+electric_sys.build_system()
+
+## Getting system
+sys, b = electric_sys.get_system(sparse_rhs=True)
+sol = spsolve(sys.tocsr(), b.todense())
+
+# Target solution (artificially made by setting voltage_source_data = np.array([5], dtype=complex))
+sol_target = np.array([
+    -1.66666667 + 1.66666667j,
+    -0.83333333 + 1.66666667j,
+    0.83333333 - 1.66666667j,
+    2.5 - 3.33333333j,
+    0.0 + 0.0j,
+    5.0 + 0.0j,
+    4.16666667 + 1.66666667j,
+])
+
+for _ in range(3000):
+    ## Computing gradients of squared error with respect to b
+    db = 2 * spsolve(sys.tocsr().conj().T, sol - sol_target)
+    drhs = db[b.row]
+
+    ## Backpropagate gradients from drhs to voltage_source_data
+    gradients = electric_sys.backpropagate_gradients(drhs=drhs)
+    ## Performing gradient descent on voltage_source_data
+    electric_sys.voltage_source_data = (
+        electric_sys.voltage_source_data - 0.01 * gradients.voltage_source_data
+    )
+
+    ## After updating voltage_source_data, rebuild the system to update sys and b
+    electric_sys.build_system()
+    sys, b = electric_sys.get_system(sparse_rhs=True)
+    sol = spsolve(sys.tocsr(), b.todense())
+
+## Checking whether we converged to the right solution
+np.testing.assert_allclose(electric_sys.voltage_source_data, np.array([5], dtype=complex))
+```
+
+For additional backpropagation examples, see `tests/test_gradients.py`.
diff --git a/docs/components/temporal-system-builder.md b/docs/components/temporal-system-builder.md
index 2193f79..c08fb6b 100644
--- a/docs/components/temporal-system-builder.md
+++ b/docs/components/temporal-system-builder.md
@@ -95,3 +95,104 @@ plt.savefig("intensities_res.png")
 This outputs the following graph that displays the intensity passing through the resistances:
 
 ![Intensity through resistances](../img/intensities_res.png)
+
+## Gradient Backpropagation
+
+`TemporalSystemBuilder` can backpropagate gradients from `S_init`, `S1`, `S2`, and `rhs` to component parameters.
+
+This makes it possible to optimize circuit parameters with gradient descent in loops where the solver output is compared to a target.
+
+### Example: Optimize Capacitor Values
+
+In this example, we optimize `capa_data` so the solution at `t=0.8` matches a target response.
+
+```python
+import numpy as np
+from scipy.sparse.linalg import spsolve
+from ElecSolver import TemporalSystemBuilder
+
+## Simple tetrahedron
+res_coords = np.array([[0, 2], [1, 3]], dtype=int)
+res_data = np.array([1, 1], dtype=float)
+
+coil_coords = np.array([[1, 0], [2, 3]], dtype=int)
+coil_data = np.array([1, 1], dtype=float)
+
+capa_coords = np.array([[1, 2], [3, 0]], dtype=int)
+capa_data = np.array([1, 1], dtype=float)
+
+## mutuals
+mutuals_coords = np.array([[], []], dtype=int)
+mutuals_data = np.array([], dtype=float)
+
+res_mutuals_coords = np.array([[], []], dtype=int)
+res_mutuals_data = np.array([], dtype=float)
+
+elec_sys = TemporalSystemBuilder(
+    coil_coords,
+    coil_data,
+    res_coords,
+    res_data,
+    capa_coords,
+    capa_data,
+    mutuals_coords,
+    mutuals_data,
+    res_mutuals_coords,
+    res_mutuals_data,
+)
+elec_sys.add_current_source(10, 1, 0)
+elec_sys.set_ground(0)
+elec_sys.build_system()
+
+## Getting initial condition system
+S_i, rhs = elec_sys.get_init_system(sparse_rhs=True)
+## Getting temporal systems
+S1, S2, rhs = elec_sys.get_system(sparse_rhs=True)
+## Solving initial condition system
+sol_init = spsolve(S_i.tocsr(), rhs.todense())
+
+## Time iteration with euler implicit scheme for 1 timestep
+dt = 0.8
+B = rhs * dt + S2 @ sol_init
+A = S2 + dt * S1
+sol = spsolve(A, B)
+
+# Target solution (artificially made by setting capa_data = np.array([0.1, 1], dtype=float))
+sol_target = np.array([
+    3.24786325,
+    -1.1965812,
+    -6.16809117,
+    0.61253561,
+    0.58404558,
+    -2.63532764,
+    0.0,
+    6.16809117,
+    2.10826211,
+    1.4957265,
+])
+
+for _ in range(1000):
+    ## computing gradients
+    dB = 2 * spsolve(A.T, sol - sol_target)
+    ## chain rule for gradients of capa_data (S2 appears twice in the computation graph)
+    dS2 = -(dB[S2.row] * sol[S2.col])
+    dS2 += dB[S2.row] * sol_init[S2.col]
+
+    ## Backpropagate gradients from dS2 to capa_data
+    gradients = elec_sys.backpropagate_gradients(dS2=dS2)
+    ## Change the values of capa_data using gradient descent
+    elec_sys.capa_data = elec_sys.capa_data - 0.01 * gradients.capa_data
+
+    ## After updating capa_data, rebuild the system to update S1, S2 and rhs
+    elec_sys.build_system()
+    S1, S2, rhs = elec_sys.get_system(sparse_rhs=True)
+    ## Recompute solution with euler implicit scheme
+    B = rhs * dt + S2 @ sol_init
+    A = S2 + dt * S1
+    sol = spsolve(A, B)
+
+## Checking whether we converged to the right solution
+np.testing.assert_allclose(elec_sys.capa_data, np.array([0.1, 1], dtype=float))
+```
+
+For additional backpropagation examples, see `tests/test_gradients.py`.
diff --git a/docs/index.md b/docs/index.md
index d65fb25..a28962e 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -10,9 +10,16 @@ This repository is **not** a general-purpose electrical system solver. Instead,
 - The graph-based description of an electric network
 - The corresponding sparse linear system to solve
 
-Its main goal is to provide a friendly Python interface for simulating analog electric systems. While suitable for small circuit simulations, its strength lies in its scalability: it is able to build linear systems with millions of nodes and components.
+In a very simple way, ElecSolver takes as an input the Resistances $R$, Capacitances $C$, Inductances $L$, Mutuals $M$, Current sources $I$  and Voltage sources $V$ along with the connectivity graph $G$ of the system and outputs the linear system: matrix $S$ and vector $b$ to solve in order to get the solution of the electric problem.
 
-!!! warning
+$$f_{\text{ElecSolver}}(R,C,L,M,I,V,G) = (S,b)$$
+
+It is then up to the user to choose the solver they want for solving the system $Sx=b$ and to manage the time iterations if needed for temporal problems.
+
+The main goal of ElecSolver is to provide a friendly Python interface for simulating and optimizing analog electric systems. While suitable for small circuit simulations, its strength lies in its scalability: it is able to build linear systems with millions of nodes and components.
+
+
+!!! tip
 
     ElecSolver has been designed with the following specifications in mind:
 
@@ -20,8 +27,9 @@ Its main goal is to provide a friendly Python interface for simulating analog el
     - Handle natively inductive mutuals and resistive mutuals.
     - Handle as many coupled electric systems as needed.
     - Deal with lonely nodes and lonely edges in the electric graph when the problem is still well posed.
+    - Allow backpropagation of gradients through the system for optimization purposes.
 
-!!! note
+!!! warning
 
     Non-linear components are not supported. You must manage event detection and system updates yourself.
 
@@ -47,12 +55,12 @@ conda install python-mumps
 
 ## Components
 
-The documentation is organized around the same component split as the README:
+The documentation presents the two main components of ElecSolver:
 
 - [FrequencySystemBuilder](components/frequency-system-builder.md)
 - [TemporalSystemBuilder](components/temporal-system-builder.md)
 
-Then follow the same supporting sections:
+Then presents extra features:
 
 - [Solver suggestions](solver-suggestions.md)
 - [Extra uses: Hydraulic or Thermal system modeling](extra-uses-hydraulic-or-thermal-system-modeling.md)
diff --git a/src/ElecSolver/FrequencySystemBuilder.py b/src/ElecSolver/FrequencySystemBuilder.py
index 85829f2..9772023 100644
--- a/src/ElecSolver/FrequencySystemBuilder.py
+++ b/src/ElecSolver/FrequencySystemBuilder.py
@@ -1,12 +1,12 @@
 import numpy as np
 import networkx as nx
 from scipy.sparse import coo_matrix,coo_array
-from .utils import SolutionFrequency
+from .utils import SolutionFrequency,GradientsParametersFrequency
 import warnings
 
 
 class FrequencySystemBuilder():
-    def __init__(self,impedence_coords,impedence_data,mutuals_coords,mutuals_data):
+    def __init__(self,impedence_coords,impedence_data,mutual_coords,mutual_data):
         """FrequencySystemBuilder class for building an electrical sparse system
         that can be solved by any sparse solver
         it supports all forms of complex making this class fit for non linear impedences
@@ -17,20 +17,20 @@ def __init__(self,impedence_coords,impedence_data,mutuals_coords,mutuals_data):
             impedence coordinates
         impedence_data : np.array of complex, shape = (N,)
             impedence value between points impedence_coords[:,i]
-        mutuals_coords : np.array of ints, shape = (2,M)
+        mutual_coords : np.array of ints, shape = (2,M)
             indexes of the impedence within impedence data which have a mutual
-        mutuals_data : np.array of complex, shape=(M,)
-            mutual value between impedences impedence_data[mutuals_coords[0,i]] impedence_data[mutuals_coords[1,i]]
+        mutual_data : np.array of complex, shape=(M,)
+            mutual value between impedences impedence_data[mutual_coords[0,i]] impedence_data[mutual_coords[1,i]]
             The mutual follows the order given in impedence coords
         """
         self.impedence_coords = impedence_coords
         self.impedence_data = impedence_data
-        self.mutuals_coords = mutuals_coords
-        self.mutuals_data = mutuals_data
-        self.current_sources_coords=np.zeros((2,0),dtype=int)
-        self.current_sources_data=np.array([],dtype=int)
-        self.voltage_sources_coords=np.zeros((2,0),dtype=int)
-        self.voltage_sources_data=np.array([],dtype=int)
+        self.mutual_coords = mutual_coords
+        self.mutual_data = mutual_data
+        self.current_source_coords=np.zeros((2,0),dtype=int)
+        self.current_source_data=np.array([],dtype=int)
+        self.voltage_source_coords=np.zeros((2,0),dtype=int)
+        self.voltage_source_data=np.array([],dtype=int)
         self.source_count = 0
         ## initializing second member as empty
         self.rhs = (np.array([]),(np.array([],dtype=int),))
@@ -38,7 +38,7 @@ def __init__(self,impedence_coords,impedence_data,mutuals_coords,mutuals_data):
         self.analysed=False
 
     def graph_analysis(self):
-        self.all_coords = np.concatenate((self.impedence_coords,self.voltage_sources_coords),axis=1)
+        self.all_coords = np.concatenate((self.impedence_coords,self.voltage_source_coords),axis=1)
         all_points = np.unique(self.all_coords)
         if all_points.shape != np.max(self.all_coords)+1:
             warnings.warn("Warning: There is one or multiple lonely nodes please clean your impedence graph")
@@ -46,7 +46,7 @@ def graph_analysis(self):
         if self.analysed:
             warnings.warn("Warning: analysis was already performed: grounds will be reasigned")
 
-        self.all_impedences = np.concatenate([self.impedence_data,self.voltage_sources_data],axis=0)
+        self.all_impedences = np.concatenate([self.impedence_data,self.voltage_source_data],axis=0)
 
         # actual number of node in the system
         self.size = np.max(self.all_coords)+1
@@ -168,9 +168,9 @@ def build_system(self):
         sign = np.sign(self.impedence_coords[0]-self.impedence_coords[1])
 
 
-        i_s_additionnal = self.offset_i + np.concatenate((self.mutuals_coords[0],self.mutuals_coords[1]),axis=0)
-        j_s_additionnal = np.concatenate((self.mutuals_coords[1],self.mutuals_coords[0]),axis=0)
-        data_additionnal = np.tile(self.mutuals_data*sign[self.mutuals_coords[0]]*sign[self.mutuals_coords[1]],(2,))
+        i_s_additionnal = self.offset_i + np.concatenate((self.mutual_coords[0],self.mutual_coords[1]),axis=0)
+        j_s_additionnal = np.concatenate((self.mutual_coords[1],self.mutual_coords[0]),axis=0)
+        data_additionnal = np.tile(self.mutual_data*sign[self.mutual_coords[0]]*sign[self.mutual_coords[1]],(2,))
 
 
         ## ground equations (1 per subsystem)
@@ -179,15 +179,25 @@ def build_system(self):
         data_ground = np.ones(len(self.affected_potentials))
 
         ## voltage sources equations
-        i_s_sources = np.tile(np.arange(0,self.voltage_sources_data.shape[0]),2)+self.number_intensities+self.size - self.source_count
-        j_s_sources = np.concatenate((self.offset_j + self.voltage_sources_coords[0],self.offset_j+self.voltage_sources_coords[1]),axis=0)
-        data_source = np.concatenate((np.ones_like(self.voltage_sources_data),-np.ones_like(self.voltage_sources_data)))
-
+        i_s_sources = np.tile(np.arange(0,self.voltage_source_data.shape[0]),2)+self.number_intensities+self.size - self.source_count
+        j_s_sources = np.concatenate((self.offset_j + self.voltage_source_coords[0],self.offset_j+self.voltage_source_coords[1]),axis=0)
+        data_source = np.concatenate((np.ones_like(self.voltage_source_data),-np.ones_like(self.voltage_source_data)))
 
+        ## Building the sparse system
         i_s = np.concatenate((i_s_nodes,i_s_edges,i_s_additionnal,i_s_ground,i_s_sources),axis=0)
         j_s = np.concatenate((j_s_nodes,j_s_edges,j_s_additionnal,j_s_ground,j_s_sources),axis=0)
         data = np.concatenate((data_nodes,data_edges,data_additionnal,data_ground,data_source),axis=0)
 
+        ##informations for gradients
+        ## building reverse S system for gradients
+        self.S_reverse_impedence = data_nodes.shape[0]+self.impedence_data.shape[0]*2 + np.arange(self.impedence_data.shape[0])
+        self.S_reverse_impedence_sign = np.ones_like(self.impedence_data)
+        # mutuals are used twice in the system but only contribute once to the gradient
+        self.S_reverse_mutual_1 = data_nodes.shape[0]+data_edges.shape[0] + np.arange(self.mutual_data.shape[0])
+        self.S_reverse_mutual_2 = data_nodes.shape[0]+data_edges.shape[0] + np.arange(self.mutual_data.shape[0],2*self.mutual_data.shape[0])
+        self.S_reverse_mutual_sign_1 = sign[self.mutual_coords[0]]*sign[self.mutual_coords[1]]
+        self.S_reverse_mutual_sign_2 = sign[self.mutual_coords[0]]*sign[self.mutual_coords[1]]
+
         self.system = (data,(i_s.astype(int),j_s.astype(int)))
         return self.system
 
@@ -211,9 +221,9 @@ def build_second_member(self,check=True):
         """
         ## Testing current
         if check:
-            for i in range(self.current_sources_coords.shape[1]):
-                input_node=self.current_sources_coords[0,i]
-                output_node=self.current_sources_coords[1,i]
+            for i in range(self.current_source_coords.shape[1]):
+                input_node=self.current_source_coords[0,i]
+                output_node=self.current_source_coords[1,i]
                 if self.number_of_subsystems>=2:
                     valid = False
                     for system in self.list_of_subgraphs:
@@ -226,10 +236,10 @@ def build_second_member(self,check=True):
 
 
         ## Building current injection
-        in_current_nodes = self.current_sources_coords[0]
-        in_current_data = - self.current_sources_data
-        out_current_nodes = self.current_sources_coords[1]
-        out_current_data = self.current_sources_data
+        in_current_nodes = self.current_source_coords[0]
+        in_current_data = - self.current_source_data
+        out_current_nodes = self.current_source_coords[1]
+        out_current_data = self.current_source_data
         current_nodes = np.concatenate((in_current_nodes,out_current_nodes),axis=0)
         current_data = np.concatenate((in_current_data,out_current_data),axis=0)
 
@@ -240,13 +250,13 @@ def build_second_member(self,check=True):
         current_nodes = current_nodes + self.rescaler[current_nodes]
 
         ## Building voltage rhs
-        voltage_nodes = np.arange(0,self.voltage_sources_data.shape[0])+self.number_intensities+self.size - self.source_count
-        voltage_data = self.voltage_sources_data
+        voltage_nodes = np.arange(0,self.voltage_source_data.shape[0])+self.number_intensities+self.size - self.source_count
+        voltage_data = self.voltage_source_data
 
         self.rhs=(np.concatenate((current_data,voltage_data),axis=0),(np.concatenate((current_nodes,voltage_nodes),axis=0),))
         return self.rhs
 
-    def get_system(self,sparse_rhs=False):
+    def get_system(self,sparse_rhs=True):
         """Function to get the system
         Parameters:
         -------
@@ -285,8 +295,8 @@ def add_current_source(self,intensity,input_node,output_node):
         output_node : int
             which node for current retrieval
         """
-        self.current_sources_coords = np.append(self.current_sources_coords,np.array([[input_node],[output_node]]),axis=1)
-        self.current_sources_data = np.append(self.current_sources_data,np.array([intensity]))
+        self.current_source_coords = np.append(self.current_source_coords,np.array([[input_node],[output_node]]),axis=1)
+        self.current_source_data = np.append(self.current_source_data,np.array([intensity]))
 
     def add_voltage_source(self,voltage,input_node,output_node):
         """Adding a voltage source to the system. This adds one equation and one degree of freedom in the system (the source intensity)
@@ -303,21 +313,94 @@ def add_voltage_source(self,voltage,input_node,output_node):
         """
         if self.analysed == True:
             warnings.warn("Warning: adding a tension source when analysis is performed may result in system topology change. You may need to rerun graph_analysis if it is the case.")
-        self.voltage_sources_coords = np.append(self.voltage_sources_coords,np.array([[input_node],[output_node]]),axis=1)
-        self.voltage_sources_data = np.append(self.voltage_sources_data,np.array([voltage]))
+        self.voltage_source_coords = np.append(self.voltage_source_coords,np.array([[input_node],[output_node]]),axis=1)
+        self.voltage_source_data = np.append(self.voltage_source_data,np.array([voltage]))
         self.source_count+=1
 
+    def backpropagate_gradients(self, dS=None, drhs=None):
+        """Function to backpropagate the gradient from the system gradient to the different parameters
+           It needs to called after the build_system function to be able to propagate the gradient on the parameters
+           For efficient backpropagation the gradients provided to this function should only be the same data arrays than S and rhs (and not the whole coo_matrix or the whole rhs vector) to avoid unnecessary computations on zero values.
+           The function will return the gradients on the parameters in the same order than they were given in the constructor of the class.
+           The returned gradients are numpy arrays of the same shape than the data arrays of the parameters.
+           For example if impedence_data was given as a parameter then the first returned gradient will be an array of the same shape than impedence_data containing the gradient on each coil value.
+           The user can provide None for the gradient that are not available or not useful to the user, in case of multiple parameters contributing to the same value in the system the function will only return the sum of the gradients on this value.
+        Parameters:
+        dS: Optional[np.array]
+            gradient on the data array of the system matrix
+        drhs: Optional[np.array]
+            gradient on the data array of rhs
+        """
+        ## Initializing gradients on parameters
+        grads_impedence = np.zeros_like(self.impedence_data,dtype=complex)
+        grads_voltage_sources = np.zeros_like(self.voltage_source_data,dtype=complex)
+        grads_mutual = np.zeros_like(self.mutual_data,dtype=complex)
+        grads_current_sources = np.zeros_like(self.current_source_data,dtype=complex)
+
+        if dS is not None:
+            ## resistance contribution to S1
+            grads_impedence += dS[self.S_reverse_impedence]*self.S_reverse_impedence_sign
+            ## resistive mutuals contribution to S1
+            grads_mutual += dS[self.S_reverse_mutual_1]*self.S_reverse_mutual_sign_1
+            grads_mutual += dS[self.S_reverse_mutual_2]*self.S_reverse_mutual_sign_2
+
+
+        if drhs is not None:
+            ## Need to make a specific treatment for current source since some equations are removed from the system
+            in_current_nodes = self.current_source_coords[0]
+            out_current_nodes = self.current_source_coords[1]
+            in_current_sign = -np.ones_like(in_current_nodes)
+            out_current_sign = np.ones_like(out_current_nodes)
+            mask_eq_in = ~np.isin(in_current_nodes,self.deleted_equation_current)
+            mask_eq_out = ~np.isin(out_current_nodes,self.deleted_equation_current)
+            in_current_nodes = in_current_nodes[mask_eq_in]
+            in_current_nodes = in_current_nodes + self.rescaler[in_current_nodes]
+            in_current_sign = in_current_sign[mask_eq_in]
+            out_current_nodes = out_current_nodes[mask_eq_out]
+            out_current_nodes = out_current_nodes + self.rescaler[out_current_nodes]
+            out_current_sign = out_current_sign[mask_eq_out]
+
+            grads_current_sources[mask_eq_in] += drhs[np.arange(in_current_nodes.shape[0])]*in_current_sign
+            grads_current_sources[mask_eq_out] += drhs[in_current_nodes.shape[0] + np.arange(out_current_nodes.shape[0])]*out_current_sign
+            grads_voltage_sources += drhs[(in_current_nodes.shape[0] + out_current_nodes.shape[0]):]
+        return GradientsParametersFrequency(grads_impedence, grads_mutual, grads_voltage_sources, grads_current_sources)
+
+
 
     def build_intensity_and_voltage_from_vector(self,sol):
+        """Function to build a more readable solution from the raw solution of the system"""
         sign = np.sign(self.all_coords[1]-self.all_coords[0])
         if self.source_count!=0:
-            return SolutionFrequency(sol[:self.number_intensities-self.source_count]*sign[:self.number_intensities-self.source_count],
-                    sol[self.number_intensities:],
+            return SolutionFrequency(sol[...,:self.number_intensities-self.source_count]*sign[:self.number_intensities-self.source_count],
+                    sol[...,self.number_intensities:],
                     sol[...,self.number_intensities-self.source_count:self.number_intensities]*sign[self.number_intensities-self.source_count:self.number_intensities]
                     )
 
         else:
-            return SolutionFrequency(sol[:self.number_intensities]*sign,
-                    sol[self.number_intensities:],
+            return SolutionFrequency(sol[...,:self.number_intensities]*sign,
+                    sol[...,self.number_intensities:],
                     np.array([],dtype=float)
-                    )
\ No newline at end of file
+                    )
+    def build_vector_from_intensity_and_voltage(self,solution: SolutionFrequency):
+        """Function to build a solution vector from a SolutionFrequency named tuple, useful for building target vectors for gradients backpropagation
+
+        Parameters
+        ----------
+        solution : SolutionFrequency
+            object containing all the components of the solution
+
+        Returns
+        -------
+        sol: np.array of shape (*, self.number_intensities+self.size)
+            array containing as many solutions as wanted
+        """
+        sign = np.sign(self.all_coords[1]-self.all_coords[0])
+        if self.source_count!=0:
+            return np.concatenate((solution.intensities*sign[:self.number_intensities-self.source_count],
+                    solution.potentials,
+                    solution.intensities_sources*sign[self.number_intensities-self.source_count:self.number_intensities]
+                    ),axis=solution.intensities.ndim-1)
+        else:
+            return np.concatenate((solution.intensities*sign,
+                    solution.potentials
+                    ),axis=solution.intensities.ndim-1)
\ No newline at end of file
diff --git a/src/ElecSolver/NetlistDumper.py b/src/ElecSolver/NetlistDumper.py
index b2a6500..141685a 100644
--- a/src/ElecSolver/NetlistDumper.py
+++ b/src/ElecSolver/NetlistDumper.py
@@ -2,7 +2,7 @@
 from .NetlistGenerator import ResistorNetlist, CapacitorNetlist, InductorNetlist, MutualInductance, SubCircuitDefinition
 from .TemporalSystemBuilder import TemporalSystemBuilder
 
-class TemporalNetlistDumper(object): 
+class TemporalNetlistDumper(object):
     def __init__(self, temporal_system_builder: TemporalSystemBuilder):
         """Constructor of TemporalNetlistDumper. Allows dumping TemporalSystemBuilder as subckt compatible with spice-like tools
 
@@ -12,15 +12,15 @@ def __init__(self, temporal_system_builder: TemporalSystemBuilder):
         if not isinstance(temporal_system_builder, TemporalSystemBuilder):
             raise TypeError("temporal_system_builder must be an instance of "+ TemporalSystemBuilder.__name__)
         self.system_builder = temporal_system_builder
-        self.ports_map ={} 
-    
+        self.ports_map ={}
+
     def add_port(self, port_number: int, port_name: str):
         """Add port to subckt definition, will be exposed and accessible in spice tools
 
         Args:
             port_number (int): port number used in temporal_system_builder to be used as port
-            port_name (str): port name to be used for export 
-        """ 
+            port_name (str): port name to be used for export
+        """
         self.ports_map.update({port_number:port_name})
 
     def generate_subcircuit_file(self, subcircuit_name : str = "esolvsub", file_name : str ="esolvsub.net"):
@@ -30,7 +30,7 @@ def generate_subcircuit_file(self, subcircuit_name : str = "esolvsub", file_name
             subcircuit_name (str, optional): name of the subcircuit to be generated. Defaults to "esolvsub".
             file_name (str, optional): file name to dump the generate the netlist to. Defaults to "esolvsub.net".
         """
-        if not bool(self.ports_map):  # check if dict is empty 
+        if not bool(self.ports_map):  # check if dict is empty
             raise ValueError("ports_map is empty, cannot generate subckt")
         list_res =[]
         list_coil=[]
@@ -60,10 +60,10 @@ def generate_subcircuit_file(self, subcircuit_name : str = "esolvsub", file_name
                 else:
                     node2=str(node2+1) # +1 because 0 node is gnd in spice !
                 curr_list_out.append(current_class(str(index), node1, node2,  value))
-        
-        for index_mutual, mutual in enumerate(self.system_builder.inductive_mutuals_data):
-            l1 = self.system_builder.inductive_mutuals_coords[0][index_mutual]
-            l2 = self.system_builder.inductive_mutuals_coords[1][index_mutual]
+
+        for index_mutual, mutual in enumerate(self.system_builder.inductive_mutual_data):
+            l1 = self.system_builder.inductive_mutual_coords[0][index_mutual]
+            l2 = self.system_builder.inductive_mutual_coords[1][index_mutual]
             l1_val = self.system_builder.coil_data[l1]
             l2_val = self.system_builder.coil_data[l2]
             k = mutual/np.sqrt(l1_val*l2_val)
diff --git a/src/ElecSolver/TemporalSystemBuilder.py b/src/ElecSolver/TemporalSystemBuilder.py
index 0c428fd..a6af190 100644
--- a/src/ElecSolver/TemporalSystemBuilder.py
+++ b/src/ElecSolver/TemporalSystemBuilder.py
@@ -1,11 +1,11 @@
 import numpy as np
 import networkx as nx
 from scipy.sparse import coo_matrix, block_diag, coo_array
-from .utils import SolutionTemporal
+from .utils import SolutionTemporal,GradientsParametersTemporal
 import warnings
 
 class TemporalSystemBuilder():
-    def __init__(self,coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,inductive_mutuals_coords,inductive_mutuals_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):
         """Class for building the linear system to solve for both temporal and frequency studies
 
         Parameters
@@ -25,18 +25,18 @@ def __init__(self,coil_coords,coil_data,res_coords,res_data,capa_coords,capa_dat
             Repetitions will be considered as another capacity in //
         capa_data : np.array of shape (R,) R being the number of capacities
             Values of capacity whose coordinates are in res_coords
-        inductive_mutuals_coords : np.array of shape (M,2) M being the number of mutuals in the system
+        inductive_mutual_coords : np.array of shape (M,2) M being the number of mutuals in the system
             The coordinates are the indices of the coils in the coil_data array that have a mutual effect
             Inductive mutuals are only supported between coils
             Repetitions will be considered as another mutual
-        inductive_mutuals_data : np.array of shape (M,) M being the number of mutuals in the system
-            Values of mutuals whose coordinates are in inductive_mutuals_coords.
+        inductive_mutual_data : np.array of shape (M,) M being the number of mutuals in the system
+            Values of mutuals whose coordinates are in inductive_mutual_coords.
             The mutual follows the sign of the nodes given in the coil_coords array
         res_mutual_coords :np.array of shape (N,2) N being the number of resistive mutuals in the system
             the coordinates are the indices of the coils or resistance in the coil_data and res_data array that have a resistive mutual effect
             resitive mutuals are only supported between coils and resistances
         res_mutual_data : np.array of shape (N,) N being the number of resistive mutuals in the system
-            Values of mutuals whose coordinates are in res_mutuals_coords.
+            Values of mutuals whose coordinates are in res_mutual_coords.
             The mutual follows the sign of the nodes given in the coil_coords or res_coord array
         """
 
@@ -46,14 +46,15 @@ def __init__(self,coil_coords,coil_data,res_coords,res_data,capa_coords,capa_dat
         self.res_data=res_data
         self.capa_coords=capa_coords
         self.capa_data=capa_data
-        self.inductive_mutuals_coords=inductive_mutuals_coords
-        self.inductive_mutuals_data=inductive_mutuals_data
+        self.inductive_mutual_coords=inductive_mutual_coords
+        self.inductive_mutual_data=inductive_mutual_data
         self.res_mutual_coords=res_mutual_coords
         self.res_mutual_data=res_mutual_data
-        self.current_sources_coords=np.zeros((2,0),dtype=int)
-        self.current_sources_data=np.array([],dtype=int)
-        self.voltage_sources_coords=np.zeros((2,0),dtype=int)
-        self.voltage_sources_data=np.array([],dtype=int)
+        ## Array for containing source information (cannot be given at initialization since it needs checks)
+        self.current_source_coords=np.zeros((2,0),dtype=int)
+        self.current_source_data=np.array([],dtype=int)
+        self.voltage_source_coords=np.zeros((2,0),dtype=int)
+        self.voltage_source_data=np.array([],dtype=int)
         self.source_count = 0
 
         ## initializing second member as empty
@@ -67,10 +68,10 @@ def graph_analysis(self):
         We test it with networkx and build some internal variables that are necessary for the system building
         In case you want to change the connectivity graph of your system you need to rerun the analysis to make check wether the topology changed or not
         """
-        if self.voltage_sources_data.shape==0 and self.current_sources_data.shape==0:
+        if self.voltage_source_data.shape==0 and self.current_source_data.shape==0:
             raise RuntimeError("The system does not have any sources (voltage or current) defined before running the graph analysis. Please define all your sources before calling set_ground, build_system or graph_analysis")
 
-        self.all_coords = np.concatenate((self.coil_coords,self.res_coords,self.capa_coords,self.voltage_sources_coords),axis=1)
+        self.all_coords = np.concatenate((self.coil_coords,self.res_coords,self.capa_coords,self.voltage_source_coords),axis=1)
         all_points = np.unique(self.all_coords)
         if all_points.shape != np.max(self.all_coords)+1:
             warnings.warn("Warning: There is one or multiple lonely nodes please clean your impedence graph")
@@ -78,7 +79,7 @@ def graph_analysis(self):
         if self.analysed:
             warnings.warn("Warning: analysis was already performed: grounds will be reasigned")
 
-        self.all_impedences = np.concatenate([self.coil_data,self.res_data,self.capa_data,self.voltage_sources_data],axis=0)
+        self.all_impedences = np.concatenate([self.coil_data,self.res_data,self.capa_data,self.voltage_source_data],axis=0)
 
         # actual number of node in the system
         self.size = np.max(self.all_coords)+1
@@ -240,9 +241,9 @@ def build_system(self):
         sign = np.sign(self.all_coords[0]-self.all_coords[1])
 
         ## inductive mutuals (contribution to S2 only)
-        i_s_additionnal_S2 = self.offset_i + np.concatenate((self.inductive_mutuals_coords[0],self.inductive_mutuals_coords[1]),axis=0)
-        j_s_additionnal_S2 = np.concatenate((self.inductive_mutuals_coords[1],self.inductive_mutuals_coords[0]),axis=0)
-        data_additionnal_S2 = np.tile(self.inductive_mutuals_data*sign[self.inductive_mutuals_coords[0]]*sign[self.inductive_mutuals_coords[1]],(2,))
+        i_s_additionnal_S2 = self.offset_i + np.concatenate((self.inductive_mutual_coords[0],self.inductive_mutual_coords[1]),axis=0)
+        j_s_additionnal_S2 = np.concatenate((self.inductive_mutual_coords[1],self.inductive_mutual_coords[0]),axis=0)
+        data_additionnal_S2 = np.tile(self.inductive_mutual_data*sign[self.inductive_mutual_coords[0]]*sign[self.inductive_mutual_coords[1]],(2,))
 
         ## resistive mutuals (contribution to S1 only)
         i_s_additionnal_S1= self.offset_i + np.concatenate((self.res_mutual_coords[0],self.res_mutual_coords[1]),axis=0)
@@ -256,9 +257,9 @@ def build_system(self):
         data_ground_S1 = np.ones(len(self.affected_potentials))
 
         ## voltage sources equations
-        i_s_sources_S1 = np.tile(np.arange(0,self.voltage_sources_data.shape[0]),2)+self.number_intensities+self.size - self.source_count
-        j_s_sources_S1 = np.concatenate((self.offset_j + self.voltage_sources_coords[0],self.offset_j+self.voltage_sources_coords[1]),axis=0)
-        data_source_S1 = np.concatenate((np.ones_like(self.voltage_sources_data),-np.ones_like(self.voltage_sources_data)))
+        i_s_sources_S1 = np.tile(np.arange(0,self.voltage_source_data.shape[0]),2)+self.number_intensities+self.size - self.source_count
+        j_s_sources_S1 = np.concatenate((self.offset_j + self.voltage_source_coords[0],self.offset_j+self.voltage_source_coords[1]),axis=0)
+        data_source_S1 = np.concatenate((np.ones_like(self.voltage_source_data),-np.ones_like(self.voltage_source_data)))
 
 
         ## building S1 system
@@ -266,16 +267,44 @@ def build_system(self):
         j_s_S1 = np.concatenate((j_s_nodes_S1,j_s_edges_coil_S1,j_s_edges_res_S1,j_s_edges_capa_S1,j_s_additionnal_S1,j_s_ground_S1,j_s_sources_S1),axis=0)
         data_S1 = np.concatenate((data_nodes_S1,data_edges_coil_S1,data_edges_res_S1,data_edges_capa_S1,data_additionnal_S1,data_ground_S1,data_source_S1),axis=0)
 
+        ## building reverse S1 system for gradients
+        self.S1_reverse_res = data_nodes_S1.shape[0]+data_edges_coil_S1.shape[0] + self.res_data.shape[0]*2 + np.arange(self.res_data.shape[0])
+        self.S1_reverse_res_sign = np.ones_like(self.res_data)
+        # resistive mutuals are used twice in the system but only contribute once to the gradient
+        self.S1_reverse_mutual_res_1 = data_nodes_S1.shape[0]+data_edges_coil_S1.shape[0] + data_edges_res_S1.shape[0] + data_edges_capa_S1.shape[0] + np.arange(self.res_mutual_data.shape[0])
+        self.S1_reverse_mutual_res_2 = data_nodes_S1.shape[0]+data_edges_coil_S1.shape[0] + data_edges_res_S1.shape[0] + data_edges_capa_S1.shape[0] + np.arange(self.res_mutual_data.shape[0],self.res_mutual_data.shape[0]*2)
+        self.S1_reverse_mutual_res_sign_1 = sign[self.res_mutual_coords[0]]*sign[self.res_mutual_coords[1]]
+        self.S1_reverse_mutual_res_sign_2 = self.S1_reverse_mutual_res_sign_1
+
         ## building S2 system
         i_s_S2 = np.concatenate((i_s_edges_coil_S2,i_s_edges_capa_S2,i_s_additionnal_S2),axis=0)
         j_s_S2 = np.concatenate((j_s_edges_coil_S2,j_s_edges_capa_S2,j_s_additionnal_S2),axis=0)
         data_S2 = np.concatenate((data_edges_coil_S2,data_edges_capa_S2,data_additionnal_S2),axis=0)
 
+        ## building reverse S2 system for gradients
+        self.S2_reverse_coil = np.arange(data_edges_coil_S2.shape[0])
+        self.S2_reverse_coil_sign = np.ones_like(data_edges_coil_S2)
+        # capacities are used twice in the system but only contribute once to the gradient
+        self.S2_reverse_capa_1 = data_edges_coil_S2.shape[0] + np.arange(self.capa_data.shape[0])
+        self.S2_reverse_capa_2 = data_edges_coil_S2.shape[0] + np.arange(self.capa_data.shape[0], self.capa_data.shape[0]*2)
+        self.S2_reverse_capa_sign_1 = np.ones_like(self.capa_data)
+        self.S2_reverse_capa_sign_2 = -np.ones_like(self.capa_data)
+        # inductive mutuals are used twice in the system but only contribute once to the gradient
+        self.S2_reverse_mutual_inductive_1 = data_edges_coil_S2.shape[0] + data_edges_capa_S2.shape[0] + np.arange(self.inductive_mutual_data.shape[0])
+        self.S2_reverse_mutual_inductive_2 = data_edges_coil_S2.shape[0] + data_edges_capa_S2.shape[0] + np.arange(self.inductive_mutual_data.shape[0], self.inductive_mutual_data.shape[0]*2)
+        self.S2_reverse_mutual_inductive_sign_1 = sign[self.inductive_mutual_coords[0]]*sign[self.inductive_mutual_coords[1]]
+        self.S2_reverse_mutual_inductive_sign_2 = self.S2_reverse_mutual_inductive_sign_1
+
         ## building S_init system
         i_s_init = np.concatenate((i_s_nodes_S1,i_s_edges_coil_S2,i_s_edges_res_S1,i_s_edges_capa_S2,i_s_ground_S1,i_s_sources_S1),axis=0)
         j_s_init = np.concatenate((j_s_nodes_S1,j_s_edges_coil_S2,j_s_edges_res_S1,j_s_edges_capa_S2,j_s_ground_S1,j_s_sources_S1),axis=0)
         data_init = np.concatenate((data_nodes_S1,data_edges_coil_S2,data_edges_res_S1,data_edges_capa_S2,data_ground_S1,data_source_S1),axis=0)
 
+        ## building reverse S_init system for gradients only resistances contribute to the gradient in S_init
+        ## inductance and capacities even though present are in zero equations
+        self.S_init_reverse_res = data_nodes_S1.shape[0]+data_edges_coil_S2.shape[0] + self.res_data.shape[0]*2 + np.arange(self.res_data.shape[0])
+        self.S_init_reverse_res_sign = np.ones_like(self.res_data)
+
         self.S_init=(data_init,(i_s_init,j_s_init))
 
         self.S1 = (data_S1,(i_s_S1.astype(int),j_s_S1.astype(int)))
@@ -303,9 +332,9 @@ def build_second_member(self,check=True):
         """
         ## Testing current
         if check:
-            for i in range(self.current_sources_coords.shape[1]):
-                input_node=self.current_sources_coords[0,i]
-                output_node=self.current_sources_coords[1,i]
+            for i in range(self.current_source_coords.shape[1]):
+                input_node=self.current_source_coords[0,i]
+                output_node=self.current_source_coords[1,i]
                 if self.number_of_subsystems>=2:
                     valid = False
                     for system in self.list_of_subgraphs:
@@ -318,10 +347,10 @@ def build_second_member(self,check=True):
 
 
         ## Building current injection
-        in_current_nodes = self.current_sources_coords[0]
-        in_current_data = - self.current_sources_data
-        out_current_nodes = self.current_sources_coords[1]
-        out_current_data = self.current_sources_data
+        in_current_nodes = self.current_source_coords[0]
+        in_current_data = - self.current_source_data
+        out_current_nodes = self.current_source_coords[1]
+        out_current_data = self.current_source_data
         current_nodes = np.concatenate((in_current_nodes,out_current_nodes),axis=0)
         current_data = np.concatenate((in_current_data,out_current_data),axis=0)
 
@@ -332,8 +361,8 @@ def build_second_member(self,check=True):
         current_nodes = current_nodes + self.rescaler[current_nodes]
 
         ## Building voltage rhs
-        voltage_nodes = np.arange(0,self.voltage_sources_data.shape[0])+self.number_intensities+self.size - self.source_count
-        voltage_data = self.voltage_sources_data
+        voltage_nodes = np.arange(0,self.voltage_source_data.shape[0])+self.number_intensities+self.size - self.source_count
+        voltage_data = self.voltage_source_data
 
         self.rhs=(np.concatenate((current_data,voltage_data),axis=0),(np.concatenate((current_nodes,voltage_nodes),axis=0),))
         return self.rhs
@@ -350,8 +379,8 @@ def add_current_source(self,intensity,input_node,output_node):
         output_node : int
             which node for current retrieval
         """
-        self.current_sources_coords = np.append(self.current_sources_coords,np.array([[input_node],[output_node]]),axis=1)
-        self.current_sources_data = np.append(self.current_sources_data,np.array([intensity]))
+        self.current_source_coords = np.append(self.current_source_coords,np.array([[input_node],[output_node]]),axis=1)
+        self.current_source_data = np.append(self.current_source_data,np.array([intensity]))
 
     def add_voltage_source(self,voltage,input_node,output_node):
         """Adding a voltage source to the system. This adds one equation and one degree of freedom in the system (the source intensity)
@@ -366,13 +395,13 @@ def add_voltage_source(self,voltage,input_node,output_node):
             node from where the voltage is enforced
 
         """
-        if self.analysed == True:
+        if self.analysed:
             warnings.warn("Warning: adding a tension source when analysis is performed may result in system topology change. You may need to rerun graph_analysis if it is the case.")
-        self.voltage_sources_coords = np.append(self.voltage_sources_coords,np.array([[input_node],[output_node]]),axis=1)
-        self.voltage_sources_data = np.append(self.voltage_sources_data,np.array([voltage]))
+        self.voltage_source_coords = np.append(self.voltage_source_coords,np.array([[input_node],[output_node]]),axis=1)
+        self.voltage_source_data = np.append(self.voltage_source_data,np.array([voltage]))
         self.source_count+=1
 
-    def get_init_system(self,sparse_rhs=False):
+    def get_init_system(self,sparse_rhs=True):
         """Function to determine the initial state of the system as coo_matrix
         Parameters
         ----------
@@ -398,7 +427,7 @@ def get_init_system(self,sparse_rhs=False):
             np.add.at(rhs, nodes, data_rhs)
         return sys,rhs
 
-    def get_system(self,sparse_rhs=False):
+    def get_system(self,sparse_rhs=True):
         """Function to get the whole system to solve as coo_matrix
         Parameters
         ----------
@@ -426,6 +455,76 @@ def get_system(self,sparse_rhs=False):
             np.add.at(rhs, nodes, data_rhs)
         return sys1,sys2,rhs
 
+    def backpropagate_gradients(self, dS1=None, dS2=None, dS_init=None, drhs=None):
+        """Function to backpropagate the gradient from the system gradient to the different parameters
+           It needs to called after the build_system function to be able to propagate the gradient on the parameters
+           For efficient backpropagation the gradients provided to this function should only be the same data arrays than S1, S2, S_init and rhs (and not the whole coo_matrix or the whole rhs vector) to avoid unnecessary computations on zero values.
+           The function will return the gradients on the parameters in the same order than they were given in the constructor of the class.
+           The returned gradients are numpy arrays of the same shape than the data arrays of the parameters.
+           For example if coil_data was given as a parameter then the first returned gradient will be an array of the same shape than coil_data containing the gradient on each coil value.
+           The user can provide None for the gradient that are not available or not useful to the user, in case of multiple parameters contributing to the same value in the system the function will only return the sum of the gradients on this value.
+           For example if a resistance is present in both S1 and S_init the returned gradient on this resistance will be the sum of the gradient coming from S1 and S_init.
+        Parameters:
+        dS1: Optional[np.array]
+            gradient on the data array of S1
+        dS2: Optional[np.array]
+            gradient on the data array of S2
+        dS_init: Optional[np.array]
+            gradient on the data array of S_init
+        drhs: Optional[np.array]
+            gradient on the data array of rhs
+        """
+        ## Initializing gradients on parameters
+        grads_coils = np.zeros_like(self.coil_data,dtype=float)
+        grads_res = np.zeros_like(self.res_data,dtype=float)
+        grads_capa = np.zeros_like(self.capa_data,dtype=float)
+        grads_voltage_sources = np.zeros_like(self.voltage_source_data,dtype=float)
+        grads_inductive_mutuals = np.zeros_like(self.inductive_mutual_data,dtype=float)
+        grads_res_mutuals = np.zeros_like(self.res_mutual_data,dtype=float)
+        grads_current_sources = np.zeros_like(self.current_source_data,dtype=float)
+
+        if dS1 is not None:
+            ## resistance contribution to S1
+            grads_res += dS1[self.S1_reverse_res]*self.S1_reverse_res_sign
+            ## resistive mutuals contribution to S1
+            grads_res_mutuals += dS1[self.S1_reverse_mutual_res_1]*self.S1_reverse_mutual_res_sign_1
+            grads_res_mutuals += dS1[self.S1_reverse_mutual_res_2]*self.S1_reverse_mutual_res_sign_2
+
+
+        if dS2 is not None:
+            ## coils contribution to S2
+            grads_coils += dS2[self.S2_reverse_coil]*self.S2_reverse_coil_sign
+            ## capacities contribution to S2
+            grads_capa += dS2[self.S2_reverse_capa_1]*self.S2_reverse_capa_sign_1
+            grads_capa += dS2[self.S2_reverse_capa_2]*self.S2_reverse_capa_sign_2
+            ## inductive mutuals contribution to S2
+            grads_inductive_mutuals += dS2[self.S2_reverse_mutual_inductive_1]*self.S2_reverse_mutual_inductive_sign_1
+            grads_inductive_mutuals += dS2[self.S2_reverse_mutual_inductive_2]*self.S2_reverse_mutual_inductive_sign_2
+
+        if dS_init is not None:
+            grads_res += dS_init[self.S_init_reverse_res]*self.S_init_reverse_res_sign
+
+        if drhs is not None:
+            ## Need to make a specific treatment for current source since some equations are removed from the system
+            in_current_nodes = self.current_source_coords[0]
+            out_current_nodes = self.current_source_coords[1]
+            in_current_sign = -np.ones_like(in_current_nodes)
+            out_current_sign = np.ones_like(out_current_nodes)
+            mask_eq_in = ~np.isin(in_current_nodes,self.deleted_equation_current)
+            mask_eq_out = ~np.isin(out_current_nodes,self.deleted_equation_current)
+            in_current_nodes = in_current_nodes[mask_eq_in]
+            in_current_nodes = in_current_nodes + self.rescaler[in_current_nodes]
+            in_current_sign = in_current_sign[mask_eq_in]
+            out_current_nodes = out_current_nodes[mask_eq_out]
+            out_current_nodes = out_current_nodes + self.rescaler[out_current_nodes]
+            out_current_sign = out_current_sign[mask_eq_out]
+
+            grads_current_sources[mask_eq_in] += drhs[np.arange(in_current_nodes.shape[0])]*in_current_sign
+            grads_current_sources[mask_eq_out] += drhs[in_current_nodes.shape[0] + np.arange(out_current_nodes.shape[0])]*out_current_sign
+            grads_voltage_sources += drhs[(in_current_nodes.shape[0] + out_current_nodes.shape[0]):]
+        return GradientsParametersTemporal(grads_coils, grads_res, grads_capa, grads_inductive_mutuals, grads_res_mutuals, grads_voltage_sources, grads_current_sources)
+
+
     def get_frequency_system(self,omega,sparse_rhs=False):
         """Function to get the complex matrix for frequency studies
 
@@ -488,3 +587,34 @@ def build_intensity_and_voltage_from_vector(self,sol):
                 np.array([],dtype=float)
             )
 
+    def build_vector_from_intensity_and_voltage(self,solution:SolutionTemporal):
+        """Utility function to build the solution vector from the different components of the solution. Useful for building a target vector for gradient backpropagation
+
+        Parameters
+        ----------
+        solution : SolutionTemporal
+            object containing all the components of the solution
+
+        Returns
+        -------
+        sol: np.array of shape (*, self.number_intensities+self.size)
+            array containing as many solutions as wanted
+        """
+        sign = np.sign(self.all_coords[1]-self.all_coords[0])
+        offset_coil = self.coil_data.shape[0]
+        offset_res = self.res_data.shape[0]
+        offset_capa = self.capa_data.shape[0]
+        if self.source_count!=0:
+            return np.concatenate((solution.coil_intensities*sign[:offset_coil],
+                solution.res_intensities*sign[offset_coil:offset_coil+offset_res],
+                solution.capa_intensities*sign[offset_coil+offset_res:offset_coil+offset_res+offset_capa],
+                solution.voltages,
+                solution.source_intensities*sign[offset_coil+offset_res+offset_capa:offset_coil+offset_res+offset_capa+self.source_count]
+            ),axis=solution.voltages.ndim-1)
+        else:
+            return np.concatenate((solution.coil_intensities*sign[:offset_coil],
+                solution.res_intensities*sign[offset_coil:offset_coil+offset_res],
+                solution.capa_intensities*sign[offset_coil+offset_res:offset_coil+offset_res+offset_capa],
+                solution.voltages
+            ),axis=solution.voltages.ndim-1)
+
diff --git a/src/ElecSolver/utils.py b/src/ElecSolver/utils.py
index 62398ae..0d381e1 100644
--- a/src/ElecSolver/utils.py
+++ b/src/ElecSolver/utils.py
@@ -5,6 +5,9 @@
 SolutionFrequency = namedtuple("ComplexSolution",["intensities","potentials","intensities_sources"])
 SolutionTemporal = namedtuple("RealSolution",["intensities_coil","intensities_res","intensities_capa","potentials","intensities_sources"])
 
+GradientsParametersFrequency= namedtuple("GradientsParametersFrequency",["impedence_data","mutual_data","voltage_source_data","current_source_data"])
+GradientsParametersTemporal = namedtuple("GradientsParametersTemporal",["coil_data","res_data","capa_data","inductive_mutual_data","res_mutual_data","voltage_source_data","current_source_data"])
+
 def parallel_sum(*impedences):
     """Function to compute the graph of impedences resulting from // graphs
     works for any number of impedence graphs
diff --git a/tests/test_gradients.py b/tests/test_gradients.py
new file mode 100644
index 0000000..8ae092e
--- /dev/null
+++ b/tests/test_gradients.py
@@ -0,0 +1,421 @@
+import numpy as np
+from scipy.sparse.linalg import spsolve
+from ElecSolver import TemporalSystemBuilder, FrequencySystemBuilder
+
+
+def test_backpropagation():
+    ## Simple tetrahedron
+    res_coords  = np.array([[0,2],[1,3]],dtype=int)
+    res_data = np.array([1,1],dtype=float)
+
+    coil_coords  = np.array([[1,0],[2,3]],dtype=int)
+    coil_data = np.array([1,1],dtype=float)
+
+    capa_coords = np.array([[1,2],[3,0]],dtype=int)
+    capa_data = np.array([1,1],dtype=float)
+
+    ## total impedance
+    mutuals_coords=np.array([[0],[1]],dtype=int)
+    mutuals_data = np.array([1],dtype=float)
+
+
+    res_mutuals_coords=np.array([[0],[1]],dtype=int)
+    res_mutuals_data = np.array([1],dtype=float)
+
+    elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+    elec_sys.add_current_source(10,1,0)
+    elec_sys.set_ground(0)
+    elec_sys.build_system()
+
+    S_i,b = elec_sys.get_init_system(sparse_rhs=True)
+    S1,S2,rhs = elec_sys.get_system(sparse_rhs=True)
+
+    dS1 = np.array([3**i for i in range(len(S1.data))])
+    dS2 = np.array([3**i for i in range(len(S2.data))])
+    dS_i = np.array([3**i for i in range(len(S_i.data))])
+    drhs = np.array([3**i for i in range(len(rhs.data))])
+    elec_sys.backpropagate_gradients(dS1=dS1)
+    elec_sys.backpropagate_gradients(dS2=dS2)
+    elec_sys.backpropagate_gradients(dS_init=dS_i)
+    elec_sys.backpropagate_gradients(drhs=drhs)
+    dS1 = np.ones_like(S1.data)*0.1
+    dS2 = np.ones_like(S2.data)*0.001
+    dS_i = np.ones_like(S_i.data)*0.01
+    drhs = np.ones_like(rhs.data,dtype=float)
+    elec_sys.backpropagate_gradients(dS1=dS1,dS2=dS2,dS_init=dS_i,drhs=drhs)
+
+
+def test_optim_res():
+    ## Simple tetrahedron
+    res_coords  = np.array([[0,2],[1,3]],dtype=int)
+    res_data = np.array([1,1],dtype=float)
+
+    coil_coords  = np.array([[1,0],[2,3]],dtype=int)
+    coil_data = np.array([1,1],dtype=float)
+
+    capa_coords = np.array([[1,2],[3,0]],dtype=int)
+    capa_data = np.array([1,1],dtype=float)
+
+    ## total impedance
+    mutuals_coords=np.array([[],[]],dtype=int)
+    mutuals_data = np.array([],dtype=float)
+
+
+    res_mutuals_coords=np.array([[],[]],dtype=int)
+    res_mutuals_data = np.array([],dtype=float)
+
+    elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+    elec_sys.add_current_source(10,1,0)
+    elec_sys.set_ground(0)
+    elec_sys.build_system()
+
+    S_i,rhs = elec_sys.get_init_system(sparse_rhs=True)
+    sol =spsolve(S_i.tocsr(),rhs.todense())
+    ## Solution when res_data = np.array([4,1],dtype=float)
+    sol_target = [ 0.,0,-2.,-8.,8.,-8.,0.,8.,0.,8.]
+    for i in range(1000):
+        ## computing gradients
+        db = 2*spsolve(S_i.tocsr().T, sol - sol_target)
+        drhs = db[rhs.row]
+        dS_i = -( db[S_i.row]*sol[S_i.col])
+        gradients = elec_sys.backpropagate_gradients(dS_init=dS_i)# ,drhs=drhs)
+        elec_sys.res_data = elec_sys.res_data - 0.01*gradients.res_data
+        elec_sys.res_data[1] = 1.
+        elec_sys.build_system()
+        S_i,rhs = elec_sys.get_init_system(sparse_rhs=True)
+        sol = spsolve(S_i.tocsr(),rhs.todense())
+    ## Checking whether we converged to the right solution
+    np.testing.assert_allclose(elec_sys.res_data, np.array([4,1],dtype=float))
+
+def test_optim_capa():
+    ## Simple tetrahedron
+    res_coords  = np.array([[0,2],[1,3]],dtype=int)
+    res_data = np.array([1,1],dtype=float)
+
+    coil_coords  = np.array([[1,0],[2,3]],dtype=int)
+    coil_data = np.array([1,1],dtype=float)
+
+    capa_coords = np.array([[1,2],[3,0]],dtype=int)
+    capa_data = np.array([1,1],dtype=float)
+
+    ## total impedance
+    mutuals_coords=np.array([[],[]],dtype=int)
+    mutuals_data = np.array([],dtype=float)
+
+
+    res_mutuals_coords=np.array([[],[]],dtype=int)
+    res_mutuals_data = np.array([],dtype=float)
+
+    elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+    elec_sys.add_current_source(10,1,0)
+    elec_sys.set_ground(0)
+    elec_sys.build_system()
+
+    S_i,rhs = elec_sys.get_init_system(sparse_rhs=True)
+    S1,S2,rhs = elec_sys.get_system(sparse_rhs=True)
+    sol_init =spsolve(S_i.tocsr(),rhs.todense())
+    dt=0.8
+    B = rhs*dt+S2@sol_init
+    A = S2+dt*S1
+    sol = spsolve(A,B)
+    ## Solution when capa_data = np.array([0.1,1],dtype=float) for first timestep
+    sol_target = [ 3.24786325, -1.1965812,  -6.16809117,  0.61253561,  0.58404558, -2.63532764, 0., 6.16809117,  2.10826211,  1.4957265 ]
+    for i in range(1000):
+        ## computing gradients
+        dB = 2*spsolve(A.T, sol - sol_target)
+        ## chain rule for gradients of capa_data (S2 appears twice in the computation graph)
+        dS2 = -( dB[S2.row]*sol[S2.col])+(dB[S2.row]*sol_init[S2.col])
+        gradients = elec_sys.backpropagate_gradients(dS2=dS2)# ,drhs=drhs)
+        elec_sys.capa_data = elec_sys.capa_data - 0.01*gradients.capa_data
+        elec_sys.capa_data[1] = 1.
+        elec_sys.build_system()
+        S1,S2,rhs = elec_sys.get_system(sparse_rhs=True)
+        B = rhs*dt+S2@sol_init
+        A = S2+dt*S1
+        sol = spsolve(A,B)
+    ## Checking whether we converged to the right solution
+    np.testing.assert_allclose(elec_sys.capa_data, np.array([0.1,1],dtype=float))
+
+def test_optim_coil():
+    ## Simple tetrahedron
+    res_coords  = np.array([[0,2],[1,3]],dtype=int)
+    res_data = np.array([1,1],dtype=float)
+
+    coil_coords  = np.array([[1,0],[2,3]],dtype=int)
+    coil_data = np.array([1,1],dtype=float)
+
+    capa_coords = np.array([[1,2],[3,0]],dtype=int)
+    capa_data = np.array([1,1],dtype=float)
+
+    ## total impedance
+    mutuals_coords=np.array([[],[]],dtype=int)
+    mutuals_data = np.array([],dtype=float)
+
+
+    res_mutuals_coords=np.array([[],[]],dtype=int)
+    res_mutuals_data = np.array([],dtype=float)
+
+    elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+    elec_sys.add_current_source(10,1,0)
+    elec_sys.set_ground(0)
+    elec_sys.build_system()
+
+    S_i,rhs = elec_sys.get_init_system(sparse_rhs=True)
+    S1,S2,rhs = elec_sys.get_system(sparse_rhs=True)
+    sol_init =spsolve(S_i.tocsr(),rhs.todense())
+    dt=0.8
+    B = rhs*dt+S2@sol_init
+    A = S2+dt*S1
+    sol = spsolve(A,B)
+    ## Solution when coil_data = np.array([3,1],dtype=float) for first timestep
+    sol_target = np.array([ 1.02628285,-2.27784731,-5.57015714,-1.1257127,3.40356001,-2.15199555,0.,5.57015714, 1.72159644, 2.84730914])
+    for i in range(1000):
+        ## computing gradients
+        dB = 2*spsolve(A.T, sol - sol_target)
+        ## chain rule for gradients of coil_data (S2 appears twice in the computation graph)
+        dS2 = -( dB[S2.row]*sol[S2.col])+(dB[S2.row]*sol_init[S2.col])
+        gradients = elec_sys.backpropagate_gradients(dS2=dS2)# ,drhs=drhs)
+        elec_sys.coil_data = elec_sys.coil_data - 0.05*gradients.coil_data
+        elec_sys.coil_data[1] = 1.
+        elec_sys.build_system()
+        S1,S2,rhs = elec_sys.get_system(sparse_rhs=True)
+        B = rhs*dt+S2@sol_init
+        A = S2+dt*S1
+        sol = spsolve(A,B)
+    ## Checking whether we converged to the right solution
+    np.testing.assert_allclose(elec_sys.coil_data, np.array([3,1],dtype=float))
+
+def test_optim_mutual():
+    ## Simple tetrahedron
+    res_coords  = np.array([[0,2],[1,3]],dtype=int)
+    res_data = np.array([1,1],dtype=float)
+
+    coil_coords  = np.array([[1,0],[2,3]],dtype=int)
+    coil_data = np.array([1,1],dtype=float)
+
+    capa_coords = np.array([[1,2],[3,0]],dtype=int)
+    capa_data = np.array([1,1],dtype=float)
+
+    ## total impedance
+    inductive_mutuals_coords=np.array([[0],[1]],dtype=int)
+    inductive_mutuals_data = np.array([0.1],dtype=float)
+
+
+    res_mutuals_coords=np.array([[],[]],dtype=int)
+    res_mutuals_data = np.array([],dtype=float)
+
+    elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,inductive_mutuals_coords,inductive_mutuals_data,res_mutuals_coords,res_mutuals_data)
+    elec_sys.add_current_source(10,1,0)
+    elec_sys.set_ground(0)
+    elec_sys.build_system()
+
+    S_i,rhs = elec_sys.get_init_system(sparse_rhs=True)
+    S1,S2,rhs = elec_sys.get_system(sparse_rhs=True)
+    sol_init =spsolve(S_i.tocsr(),rhs.todense())
+    dt=0.8
+    B = rhs*dt+S2@sol_init
+    print(B)
+    A = S2+dt*S1
+    sol = spsolve(A,B)
+    ## Solution when mutuals_data = np.array([0.5],dtype=float) for first timestep
+    sol_target = np.array([ 3.07692308, -3.07692308, -4.14529915,  0.2991453,   2.77777778, -2.77777778, 0., 4.14529915,  2.22222222,  1.92307692])
+    for i in range(1000):
+        ## computing gradients
+        dB = 2*spsolve(A.T, sol - sol_target)
+        ## chain rule for gradients of coil_data (S2 appears twice in the computation graph)
+        dS2 = -( dB[S2.row]*sol[S2.col])+(dB[S2.row]*sol_init[S2.col])
+        gradients = elec_sys.backpropagate_gradients(dS2=dS2)# ,drhs=drhs)
+        elec_sys.inductive_mutual_data = elec_sys.inductive_mutual_data - 0.001*gradients.inductive_mutual_data
+        elec_sys.build_system()
+        S1,S2,rhs = elec_sys.get_system(sparse_rhs=True)
+        B = rhs*dt+S2@sol_init
+        A = S2+dt*S1
+        sol = spsolve(A,B)
+    ## Checking whether we converged to the right solution
+    np.testing.assert_allclose(elec_sys.inductive_mutual_data, np.array([0.5],dtype=float))
+
+def test_optim_res_mutual():
+    ## Simple tetrahedron
+    res_coords  = np.array([[0,2],[1,3]],dtype=int)
+    res_data = np.array([1,1],dtype=float)
+
+    coil_coords  = np.array([[1,0],[2,3]],dtype=int)
+    coil_data = np.array([1,1],dtype=float)
+
+    capa_coords = np.array([[1,2],[3,0]],dtype=int)
+    capa_data = np.array([1,1],dtype=float)
+
+    ## total impedance
+    inductive_mutuals_coords=np.array([[],[]],dtype=int)
+    inductive_mutuals_data = np.array([],dtype=float)
+
+
+    res_mutuals_coords=np.array([[0],[1]],dtype=int)
+    res_mutuals_data = np.array([0.1],dtype=float)
+
+    elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,inductive_mutuals_coords,inductive_mutuals_data,res_mutuals_coords,res_mutuals_data)
+    elec_sys.add_current_source(10,1,0)
+    elec_sys.set_ground(0)
+    elec_sys.build_system()
+
+    S_i,rhs = elec_sys.get_init_system(sparse_rhs=True)
+    S1,S2,rhs = elec_sys.get_system(sparse_rhs=True)
+    sol_init =spsolve(S_i.tocsr(),rhs.todense())
+    dt=0.8
+    B = rhs*dt+S2@sol_init
+    A = S2+dt*S1
+    sol = spsolve(A,B)
+    ## Solution when res mutuals_data = np.array([0.5],dtype=float) for first timestep
+    sol_target = np.array([ 2.85714286, -2.85714286, -4.36507937,  0.07936508,  2.77777778, -2.77777778, 0., 4.36507937,  2.22222222,  2.14285714])
+    for i in range(2000):
+        ## computing gradients
+        dB = 2*spsolve(A.T, sol - sol_target)
+        ## chain rule for gradients of coil_data (S1 appears once in the computation graph)
+        dS1 = -( dB[S1.row]*sol[S1.col])*dt
+        gradients = elec_sys.backpropagate_gradients(dS1=dS1)# ,drhs=drhs)
+        elec_sys.res_mutual_data = elec_sys.res_mutual_data - 0.05*gradients.res_mutual_data
+        elec_sys.build_system()
+        S1,S2,rhs = elec_sys.get_system(sparse_rhs=True)
+        B = rhs*dt+S2@sol_init
+        A = S2+dt*S1
+        sol = spsolve(A,B)
+    ## Checking whether we converged to the right solution
+    np.testing.assert_allclose(elec_sys.res_mutual_data, np.array([0.5],dtype=float))
+
+def test_optim_current_sources():
+    ## Simple tetrahedron
+    res_coords  = np.array([[0,2],[1,3]],dtype=int)
+    res_data = np.array([1,1],dtype=float)
+
+    coil_coords  = np.array([[1,0],[2,3]],dtype=int)
+    coil_data = np.array([1,1],dtype=float)
+
+    capa_coords = np.array([[1,2],[3,0]],dtype=int)
+    capa_data = np.array([1,1],dtype=float)
+
+    ## total impedance
+    mutuals_coords=np.array([[],[]],dtype=int)
+    mutuals_data = np.array([],dtype=float)
+
+
+    res_mutuals_coords=np.array([[],[]],dtype=int)
+    res_mutuals_data = np.array([],dtype=float)
+
+    elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+    elec_sys.add_current_source(10,1,0)
+    elec_sys.set_ground(0)
+    elec_sys.build_system()
+
+    S_i,rhs = elec_sys.get_init_system(sparse_rhs=True)
+    sol =spsolve(S_i.tocsr(),rhs.todense())
+    ## Solution when current_source_data = np.array([8],dtype=float)
+    sol_target = np.array([ 0.,  0., -4., -4.,  4., -4.,  0.,  4.,  0.,  4.])
+    for i in range(1000):
+        ## computing gradients (source information is in the rhs of the system)
+        db = 2*spsolve(S_i.tocsr().T, sol - sol_target)
+        drhs = db[rhs.row]
+        ## backpropagation of gradients to current_source_data
+        gradients = elec_sys.backpropagate_gradients(drhs=drhs)
+        elec_sys.current_source_data = elec_sys.current_source_data - 0.01*gradients.current_source_data
+        elec_sys.build_system()
+        S_i,rhs = elec_sys.get_init_system(sparse_rhs=True)
+        sol = spsolve(S_i.tocsr(),rhs.todense())
+    ## Checking whether we converged to the right solution
+    np.testing.assert_allclose(elec_sys.current_source_data, np.array([8],dtype=float))
+
+def test_backpropagation_frequency():
+    ## Simple circuit with one coil and one res
+    impedence_coords = np.array([[0,0],[1,2]],dtype=int)
+    impedence_data = np.array([1,1j],dtype=complex)
+
+    mutuals_coords=np.array([[0],[1]],dtype=int)
+    mutuals_data = np.array([0.1j],dtype=complex)
+
+    electric_sys = FrequencySystemBuilder(impedence_coords,impedence_data,mutuals_coords,mutuals_data)
+    electric_sys.add_current_source(intensity=10,input_node=1,output_node=0)
+    electric_sys.set_ground(0)
+    electric_sys.build_system()
+
+    S,b = electric_sys.get_system(sparse_rhs=True)
+    dS = np.array([3**i+1j for i in range(len(S.data))])
+    drhs = np.array([3**i+1j for i in range(len(b.data))])
+    electric_sys.backpropagate_gradients(dS=dS,drhs=drhs)
+
+def test_optim_mutual_complex():
+    ## sparse python res matrix
+    impedence_coords = np.array([[0,0,1],[1,2,2]],dtype=int)
+    impedence_data = np.array([1,1,1],dtype=complex)
+
+    ## mutuals
+    mutuals_coords=np.array([[0],[1]],dtype=int)
+    mutuals_data = np.array([2.j],dtype=complex)
+
+
+    electric_sys = FrequencySystemBuilder(impedence_coords,impedence_data,mutuals_coords,mutuals_data)
+    electric_sys.add_voltage_source(voltage=10,input_node=1,output_node=0)
+    # setting the ground
+    electric_sys.set_ground(0)
+    electric_sys.build_system()
+
+    ## Need to evaluate the system because it was altered when calling the second member
+    sys,b = electric_sys.get_system(sparse_rhs=True)
+    sol = spsolve(sys.tocsr(),b.todense())
+    ## Target solution when impedence_data = np.array([1+1j,1,1],dtype=complex)
+    sol_target = np.array([-2.+4.j, -1.+2.j,  1.-2.j,  3.-6.j,  0.+0.j, 10.+0.j,  9.+2.j])
+    for i in range(3000):
+        ## computing gradients
+        db = 2*spsolve(sys.tocsr().conj().T, sol - sol_target)
+        dS = -( db[sys.row]*np.conj(sol[sys.col]))
+        gradients = electric_sys.backpropagate_gradients(dS=dS)# ,drhs=drhs)
+        electric_sys.impedence_data = electric_sys.impedence_data - 0.01*gradients.impedence_data
+        electric_sys.build_system()
+        sys,b = electric_sys.get_system(sparse_rhs=True)
+        sol = spsolve(sys.tocsr(),b.todense())
+    ## Checking whether we converged to the right solution
+    np.testing.assert_allclose(electric_sys.impedence_data, np.array([1+1j,1,1],dtype=complex))
+
+def test_optim_source_complex():
+    ## sparse python res matrix
+    impedence_coords = np.array([[0,0,1],[1,2,2]],dtype=int)
+    impedence_data = np.array([1,1,1],dtype=complex)
+
+    ## mutuals
+    mutuals_coords=np.array([[0],[1]],dtype=int)
+    mutuals_data = np.array([2.j],dtype=complex)
+
+
+    electric_sys = FrequencySystemBuilder(impedence_coords,impedence_data,mutuals_coords,mutuals_data)
+    electric_sys.add_voltage_source(voltage=10,input_node=1,output_node=0)
+    # setting the ground
+    electric_sys.set_ground(0)
+    electric_sys.build_system()
+
+    ## Need to evaluate the system because it was altered when calling the second member
+    sys,b = electric_sys.get_system(sparse_rhs=True)
+    sol = spsolve(sys.tocsr(),b.todense())
+    ## Target solution when voltage_source_data = np.array([5],dtype=complex)
+    sol_target = np.array([-1.66666667+1.66666667j,-0.83333333+1.66666667j,0.83333333-1.66666667j,2.5       -3.33333333j,  0.        +0.j,          5.        +0.j, 4.16666667+1.66666667j])
+    for i in range(3000):
+        ## computing gradients
+        db = 2*spsolve(sys.tocsr().conj().T, sol - sol_target)
+        drhs = db[b.row]
+        gradients = electric_sys.backpropagate_gradients(drhs=drhs)
+        electric_sys.voltage_source_data = electric_sys.voltage_source_data - 0.01*gradients.voltage_source_data
+        electric_sys.build_system()
+        sys,b = electric_sys.get_system(sparse_rhs=True)
+        sol = spsolve(sys.tocsr(),b.todense())
+    ## Checking whether we converged to the right solution
+    np.testing.assert_allclose(electric_sys.voltage_source_data, np.array([5],dtype=complex))
+
+
+if __name__ == "__main__":
+    test_backpropagation()
+    test_optim_res()
+    test_optim_capa()
+    test_optim_coil()
+    test_optim_mutual()
+    test_optim_res_mutual()
+    test_optim_current_sources()
+    test_backpropagation_frequency()
+    test_optim_mutual_complex()
+    test_optim_source_complex()
\ No newline at end of file
diff --git a/tests/test_sysbuilder.py b/tests/test_sysbuilder.py
index 9a9e80f..614ca48 100644
--- a/tests/test_sysbuilder.py
+++ b/tests/test_sysbuilder.py
@@ -5,7 +5,6 @@
 
 
 def test_sys_general_coo_build():
-    jw = 5j
     ## sparse python res matrix
     impedence_coords = np.array([[0,0,1],[1,2,2]],dtype=int)
     impedence_data = np.array([1,1,1],dtype=complex)
@@ -67,12 +66,12 @@ def test_res_grid():
     mask = adjacency.row>adjacency.col
     impedence_coords = np.array([adjacency.row,adjacency.col],dtype=int)[:,mask]
     impedence_data = adjacency.data[mask]
-    mutual_coords = [[],[]]
-    mutual_data = []
+    mutual_coords = np.array([[],[]],dtype=int)
+    mutual_data = np.array([],dtype=complex)
 
 
 
-    electric_sys = FrequencySystemBuilder(impedence_coords=impedence_coords,impedence_data=impedence_data,mutuals_coords=mutual_coords,mutuals_data=mutual_data)
+    electric_sys = FrequencySystemBuilder(impedence_coords=impedence_coords,impedence_data=impedence_data,mutual_coords=mutual_coords,mutual_data=mutual_data)
     electric_sys.add_current_source(intensity=2,input_node=0,output_node=24)
     electric_sys.add_current_source(intensity=2,input_node=42,output_node=24)
     electric_sys.add_current_source(intensity=2,input_node=48,output_node=24)
diff --git a/tests/test_temporal_netlist_dumper.py b/tests/test_temporal_netlist_dumper.py
index 27213ff..714464e 100644
--- a/tests/test_temporal_netlist_dumper.py
+++ b/tests/test_temporal_netlist_dumper.py
@@ -1,4 +1,4 @@
-from ElecSolver import TemporalSystemBuilder 
+from ElecSolver import TemporalSystemBuilder
 from ElecSolver import TemporalNetlistDumper
 import numpy as np
 import unittest
@@ -17,12 +17,12 @@ def setUp(self):
         capa_coords = np.array([[1,2],[3,0]],dtype=int)
         capa_data = np.array([1,1],dtype=float)
         ## coupling factors
-        mutuals_coords=np.array([[0],[1]],dtype=int)
-        mutuals_data = np.array([1e-5],dtype=float)
-        res_mutuals_coords=np.array([[],[]],dtype=int)
-        res_mutuals_data = np.array([],dtype=float)
+        mutual_coords=np.array([[0],[1]],dtype=int)
+        mutual_data = np.array([1e-5],dtype=float)
+        res_mutual_coords=np.array([[],[]],dtype=int)
+        res_mutual_data = np.array([],dtype=float)
 
-        self.elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+        self.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)
 
     def tearDown(self):
         # Remove the test file after tests
@@ -36,7 +36,7 @@ def test_init_invalid_system(self):
 
     def test_init_valid_system(self):
         TemporalNetlistDumper(self.elec_sys)
-    
+
     def test_add_port(self):
         test = TemporalNetlistDumper(self.elec_sys)
         test.add_port(1,'in')
diff --git a/tests/test_temporal_system.py b/tests/test_temporal_system.py
index a0f34c6..f0fe821 100644
--- a/tests/test_temporal_system.py
+++ b/tests/test_temporal_system.py
@@ -19,14 +19,14 @@ def test_temporal():
     capa_data = np.array([1],dtype=float)
 
     ## total impedance
-    mutuals_coords=np.array([[],[]],dtype=int)
-    mutuals_data = np.array([],dtype=float)
+    mutual_coords=np.array([[],[]],dtype=int)
+    mutual_data = np.array([],dtype=float)
 
 
-    res_mutuals_coords=np.array([[],[]],dtype=int)
-    res_mutuals_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,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+    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_current_source(10,1,0)
     elec_sys.set_ground(0)
     elec_sys.build_system()
@@ -82,14 +82,14 @@ def test_temporal2():
     capa_data = np.array([2,2],dtype=float)
 
     ## total impedance
-    mutuals_coords=np.array([[],[]],dtype=int)
-    mutuals_data = np.array([],dtype=float)
+    mutual_coords=np.array([[],[]],dtype=int)
+    mutual_data = np.array([],dtype=float)
 
 
-    res_mutuals_coords=np.array([[],[]],dtype=int)
-    res_mutuals_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,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+    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_current_source(10,1,0)
     elec_sys.set_ground(0)
     elec_sys.build_system()
@@ -136,21 +136,21 @@ def test_one_shot_temporal():
     capa_data = np.array([1,1],dtype=float)
 
     ## total impedance
-    mutuals_coords=np.array([[],[]],dtype=int)
-    mutuals_data = np.array([],dtype=float)
+    mutual_coords=np.array([[],[]],dtype=int)
+    mutual_data = np.array([],dtype=float)
 
 
-    res_mutuals_coords=np.array([[],[]],dtype=int)
-    res_mutuals_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,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+    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_current_source(10,1,0)
     elec_sys.set_ground(0)
     elec_sys.build_system()
 
-    S_i,b = elec_sys.get_init_system()
+    S_i,b = elec_sys.get_init_system(sparse_rhs=False)
     sol = spsolve(S_i.tocsr(),b)
-    S1,S2,rhs = elec_sys.get_system()
+    S1,S2,rhs = elec_sys.get_system(sparse_rhs=False)
     dt=0.08
     nb_timesteps=100
     ## build big system for no for loop!
@@ -188,14 +188,14 @@ def test_voltage():
     capa_data = np.array([1,1],dtype=float)
 
     ## total impedance
-    mutuals_coords=np.array([[],[]],dtype=int)
-    mutuals_data = np.array([],dtype=float)
+    mutual_coords=np.array([[],[]],dtype=int)
+    mutual_data = np.array([],dtype=float)
 
 
-    res_mutuals_coords=np.array([[],[]],dtype=int)
-    res_mutuals_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,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+    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,1,0)
     elec_sys.set_ground(0)
     elec_sys.build_system()
@@ -226,10 +226,8 @@ def test_big_grid():
     coords_capa = np.array([[],[]],dtype=int)
     data_capa = np.array([],dtype=float)
 
-    mutual_coords = [[],[]]
-    mutual_data = []
-
-
+    mutual_coords = np.array([[],[]],dtype=int)
+    mutual_data = np.array([],dtype=float)
 
 
     electric_sys = TemporalSystemBuilder(coords_coil,data_coil,coords_res,data_res,coords_capa,data_capa,mutual_coords,mutual_data,mutual_coords,mutual_data)
@@ -306,15 +304,15 @@ def freq_simulation():
     capa_data = np.array([1,1],dtype=float)
 
     ## total impedance
-    mutuals_coords=np.array([[],[]],dtype=int)
-    mutuals_data = np.array([],dtype=float)
+    mutual_coords=np.array([[],[]],dtype=int)
+    mutual_data = np.array([],dtype=float)
 
 
-    res_mutuals_coords=np.array([[],[]],dtype=int)
-    res_mutuals_data = np.array([],dtype=float)
+    res_mutual_coords=np.array([[],[]],dtype=int)
+    res_mutual_data = np.array([],dtype=float)
 
     fvect = [10,20,100,1000]
-    elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+    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_current_source(10,1,0)
     elec_sys.set_ground(0)
     elec_sys.build_system()
@@ -341,16 +339,16 @@ def test_hydraulic():
     capa_coords = np.array([[],[]],dtype=int)
     capa_data = np.array([],dtype=float)
 
-    ## Defining empty mutuals here
-    mutuals_coords=np.array([[],[]],dtype=int)
-    mutuals_data = np.array([],dtype=float)
+    ## Defining empty mutual here
+    mutual_coords=np.array([[],[]],dtype=int)
+    mutual_data = np.array([],dtype=float)
 
 
-    res_mutuals_coords=np.array([[],[]],dtype=int)
-    res_mutuals_data = np.array([],dtype=float)
+    res_mutual_coords=np.array([[],[]],dtype=int)
+    res_mutual_data = np.array([],dtype=float)
 
     ## initializing system
-    hydraulic_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+    hydraulic_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutual_coords,mutual_data,res_mutual_coords,res_mutual_data)
     ## enforcing a pressure delta of 10 Pa
     hydraulic_sys.add_voltage_source(10,1,0)
     ## Seting ground at point 0
@@ -386,15 +384,15 @@ def test_lonely_nodes():
     capa_coords = np.array([[1],[3]],dtype=int)
     capa_data = np.array([2],dtype=float)
 
-    ## Defining empty mutuals here
-    mutuals_coords=np.array([[],[]],dtype=int)
-    mutuals_data = np.array([],dtype=float)
+    ## Defining empty mutual here
+    mutual_coords=np.array([[],[]],dtype=int)
+    mutual_data = np.array([],dtype=float)
 
 
-    res_mutuals_coords=np.array([[],[]],dtype=int)
-    res_mutuals_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,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
+    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_current_source(10,1,0)
     elec_sys.set_ground(0)
     elec_sys.build_system()
diff --git a/zensical.toml b/zensical.toml
index 76c3f45..bd52789 100644
--- a/zensical.toml
+++ b/zensical.toml
@@ -19,4 +19,24 @@ nav = [
 ]
 
 [project.theme]
-name = "zensical"
\ No newline at end of file
+name = "zensical"
+
+# Palette toggle for automatic mode
+[[project.theme.palette]]
+media = "(prefers-color-scheme)"
+toggle.icon = "lucide/sun-moon"
+toggle.name = "Switch to light mode"
+
+# Palette toggle for light mode
+[[project.theme.palette]]
+media = "(prefers-color-scheme: light)"
+scheme = "default"
+toggle.icon = "lucide/sun"
+toggle.name = "Switch to dark mode"
+
+# Palette toggle for dark mode
+[[project.theme.palette]]
+media = "(prefers-color-scheme: dark)"
+scheme = "slate"
+toggle.icon = "lucide/moon"
+toggle.name = "Switch to system preference"
\ No newline at end of file