diff --git a/docs/make.jl b/docs/make.jl
index 6938f0c..422ce18 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -16,6 +16,7 @@ makedocs(
             "Network" => "Models/network.md",
             "Generator" => "Models/gens.md",
             "Inverter" => "Models/inverters.md",
+            "Small Signal" => "Models/small.md",
         ],
     ],
 )
diff --git a/docs/src/Examples/example_OMIB.md b/docs/src/Examples/example_OMIB.md
index f6c747a..a248672 100644
--- a/docs/src/Examples/example_OMIB.md
+++ b/docs/src/Examples/example_OMIB.md
@@ -1,7 +1,7 @@
 # Tutorial: One Machine against Infinite Bus (OMIB)
 
 This tutorial will introduce you to the functionality of `LITS` for running Power System Simulations.
-Note that this tutorial is for `LITS 0.2.0`. Future versions will have dedicated functions to find an equilibrium point and a proper functions for running perturbations without coding directly the callbacks.
+Note that this tutorial is for `LITS 0.3.0`. Future versions will have dedicated functions to find an equilibrium point using a Power Flow method without relying in a guess of the initial condition to run a non-linear solver.
 
 This tutorial presents a simulation of a two-bus system, with an infinite bus (represented as a voltage source behind an impedance) at bus 1 and a classic machine on bus 2. The perturbation will be the trip of one of the two circuits (doubling its resistance and impedance) of the line that connects both buses.
 
@@ -14,13 +14,11 @@ The first step consists in initialize all packages that will be used to run the
 ```julia
 using LITS
 using PowerSystems
-using NLsolve
-using DiffEqBase
 using Sundials
 const PSY = PowerSystems
 ```
 
-`LITS` and `PowerSystems` are used to properly define the data structure, while `NLsolve`, `DiffEqBase` and `Sundials` are used to formulate and solve the problem. Finally we call use can call `PowerSystems` functions using the `PSY` abbreviation.
+`PowerSystems` is used to properly define the data structure, while `Sundials` is used to fsolve the problem defined in `LITS`. Finally we call use can call `PowerSystems` functions using the `PSY` abbreviation.
 
 ## Step 2: Data creation
 
@@ -28,60 +26,69 @@ Next we need to define the different elements required to run a simulation. To r
 
 - Vector of `PSY.Bus` elements, that define all the buses in the network.
 - Vector of `PSY.Branch` elements, that define all the branches elements (that connect two buses) in the network.
-- Vector of `DynInjection` elements, that define all the devices connected to buses that can inject (or withdraw) current, while also defining differential equations to model its dynamics.
-- Vector of `PSY.Injection` elements, that define all the devices connected to buses that can inject (or withdraw) current, without defining any differential equation.
+- Vector of `PSY.DynamicInjection` elements, that define all the devices connected to buses that can inject (or withdraw) current, while also defining differential equations to model its dynamics. These include generators and inverters.
+- Vector of `PSY.PowerLoad` elements, that define all the loads connected to buses that can withdraw current, without defining any differential equation. Note that `LITS` will convert ConstantPower loads to RLC loads for transient simulations.
+- Vector of `PSY.Source` elements, that define source components behind a reactance that can inject or withdraw current, without defining any differential equation.
 - The base of power used to define per unit values, in MVA as a `Float64` value.
 - The base frequency used in the system, in Hz as a `Float64` value.
-- (Optional) Vector of `DynBranch` elements, that can be used to model lines with differential equations.
+- (Optional) Selecting which of the `PSY.Lines` (of the `PSY.Branch` vector) elements must be modeled of `DynamicLines` elements, that can be used to model lines with differential equations.
 
 To start we will define the data structures for the network.
 
 ### Buses and Branches
 
-As mentioned earlier, we require to create a `Vector` of `PSY.Bus` to define the buses in the network. Currently, some of the parameters are not used in `LITS`, but will be used once the initialization procedure is implemented.
+As mentioned earlier, we require to create a `Vector` of `PSY.Bus` to define the buses in the network. Currently, some of the parameters are not used in `LITS`, but will be used once the initialization procedure is implemented (such as voltage limits or the requested bus voltage).
 
 ```julia
 #Define the vector of buses
-nodes_case1 = [PSY.Bus(1 , #number
-                   "Bus 1", #Name
-                   "REF" , #BusType (REF, PV, PQ)
-                   0, #Angle in radians
-                   1.05, #Voltage in pu
-                   (min=0.94, max=1.06), #Voltage limits in pu
-                   69), #Base voltage in kV
-                   PSY.Bus(2 , "Bus 2"  , "PV" , 0 , 1.0 , (min=0.94, max=1.06), 69)];
+nodes_OMIB = [
+    PSY.Bus(
+        1, #number
+        "Bus 1", #Name
+        "REF", #BusType (REF, PV, PQ)
+        0, #Angle in radians
+        1.05, #Voltage in pu
+        (min = 0.94, max = 1.06), #Voltage limits in pu
+        69, #Base voltage in kV
+    ),
+    PSY.Bus(2, "Bus 2", "PV", 0, 1.0, (min = 0.94, max = 1.06), 69),
+]
 ```
 
-Note that two buses are defined in the vector `nodes_case1`. Similarly, to define the branches (that also has some parameters that are currently not used):
+Note that two buses are defined in the vector `nodes_case1`. It is important that the bus numbers are ordered from ``1`` to ``n``, since that structure will be used to construct the vector of variables. Future versions of `LITS` will allow to relax this assumption. Similarly, to define the branches (that also has some parameters that are currently not used, such as the rate and angle limits):
 
 ```julia
 #Define the vector of branches
-branch_case1 = [PSY.Line("Line1", #name
-                     true, #available
-                     0.0, #active power flow initial condition (from-to)
-                     0.0, #reactive power flow initial condition (from-to)
-                     Arc(from=nodes_case1[1], to=nodes_case1[2]), #Connection between buses
-                     0.01, #resistance in pu
-                     0.05, #reactance in pu
-                     (from=0.0, to=0.0), #susceptance in pu
-                     18.046, #rate in MW
-                     1.04)];  #angle limits (-min and max)
+branch_OMIB = [PSY.Line(
+    "Line1", #name
+    true, #available
+    0.0, #active power flow initial condition (from-to)
+    0.0, #reactive power flow initial condition (from-to)
+    Arc(from = nodes_OMIB[1], to = nodes_OMIB[2]), #Connection between buses
+    0.01, #resistance in pu
+    0.05, #reactance in pu
+    (from = 0.0, to = 0.0), #susceptance in pu
+    18.046, #rate in MW
+    1.04, #angle limits (-min and max)
+)]
 ```
 
 Since we are interested in creating a fault that trips one of the two circuits of the line, we will create an additional `Vector` of branches with doubled impedance:
 
 ```julia
 #Define the vector of branches under the fault
-branch_case1_fault = [PSY.Line("Line1", #name
-                           true, #available
-                           0.0, #active power flow initial condition (from-to)
-                           0.0, #reactive power flow initial condition (from-to)
-                           Arc(from=nodes_case1[1], to=nodes_case1[2]), #Connection between buses
-                           0.02, #resistance in pu
-                           0.1, #reactance in pu
-                           (from=0.0, to=0.0), #susceptance in pu
-                           18.046, #rate in MW
-                           1.04)];  #angle limits (-min and max)
+branch_OMIB_fault = [PSY.Line(
+    "Line1", #name
+    true, #available
+    0.0, #active power flow initial condition (from-to)
+    0.0, #reactive power flow initial condition (from-to)
+    Arc(from = nodes_OMIB[1], to = nodes_OMIB[2]), #Connection between buses
+    0.02, #resistance in pu
+    0.1, #reactance in pu
+    (from = 0.0, to = 0.0), #susceptance in pu
+    18.046, #rate in MW
+    1.04, #angle limits (-min and max)
+)]
 ```
 
 Note that the resistance and reactance is doubled in comparison to the system without fault.
@@ -91,161 +98,215 @@ Note that the resistance and reactance is doubled in comparison to the system wi
 Secondly, we will define devices that can inject/withdraw electric current directly without defining differential equations. In this case we include a load and the voltage source that model the infinite bus.
 
 ```julia
-loads_case1 = [PowerLoad("Bus1", #name
-                         true, #available
-                         nodes_case1[2], #bus location
-                         PowerSystems.ConstantPower, #type
-                         0.3, #Active Power pu
-                         0.01, #Reactive power pu
-                         0.3, #Max Active Power pu
-                         0.01)]; #Max Reactive Power pu
-
-inf_gen_case1 = StaticSource(1, #number
-                :InfBus, #name
-                nodes_case1[1], #bus
-                1.05, #VR real part of voltage source
-                0.0, #VI imaginary part of voltage source
-                0.000001); #Xth
+loads_OMIB = [PSY.PowerLoad(
+    "LBus1", #name
+    true, #availability
+    nodes_OMIB[2], #bus
+    PSY.ConstantPower, #type
+    0.3, #P
+    0.01, #Q
+    0.3, #P_max
+    0.01, #Q_max
+)]
+
+inf_gen_OMIB = [PSY.Source(
+    "InfBus", #name
+    true, #availability
+    nodes_OMIB[1], #bus
+    1.05, #VR
+    0.0, #VI
+    0.000005, #Xth
+)]
 ```
 
 Note that loads are assumed as constant power for power flow purposes, but for dynamic simulations are converted to impedance loads assuming nominal voltage equals to 1 pu.
 
 ### Dynamic Injection devices
 
-Third, we define the `Vector` of `DynInjection` elements. In this case, we require to define a generator located in bus 2. For that purpose, we need to define its machine, shaft, automatic voltage regulator (AVR), turbine governor (TG) and power system stabilizer (PSS):
+Third, we define the `Vector` of `PSY.DynamicInjection` elements. In this case, we require to define a generator located in bus 2. For that purpose, we need to define its machine, shaft, automatic voltage regulator (AVR), turbine governor (TG) and power system stabilizer (PSS):
 
 ```julia
 ### Machine ###
-case1_machine = BaseMachine(0.0, #R
-                            0.2995, #Xd_p
-                            0.7087, #eq_p
-                            100.0)  #MVABase
+machine_OMIB = PSY.BaseMachine(
+    0.0, #R
+    0.2995, #Xd_p
+    0.7087, #eq_p
+    100.0,
+)  #MVABase
 
 ######## Shaft Data #########
 
 ### Shaft for Case 1 ###
-case1_shaft = SingleMass(3.148, #H
-                     2.0) #D
+shaft_OMIB = PSY.SingleMass(
+    3.148, #H
+    2.0, #D
+)
 
 ########  AVR Data #########
-case1_avr = AVRFixed(0.0) #Vf not applicable in Classic Machines
+avr_OMIB = AVRFixed(0.0) #Vf not applicable in Classic Machines
 
 ######## TG Data #########
 ### No TG ###
-case1234_no_tg = TGFixed(1.0) #eff
+tg_OMIB = TGFixed(1.0) #No TG: Efficiency = 1.0
 
 ######## PSS Data #########
 ### No PSS ###
-cases_no_pss = PSSFixed(0.0) #No PSS without AVR
+pss_OMIB = PSSFixed(0.0) #No PSS without AVR
 
 ### Constructing the Generator ###
-case1_gen = PSY.DynamicGenerator(1, #Number
-                      :Case1Gen,
-                      nodes_case1[2], #bus
-                      1.0, # ω_ref,
-                      1.0, #V_ref (only used in AVR)
-                      0.5, #P_ref
-                      case1_machine, #machine
-                      case1_shaft, #shaft
-                      case1_avr, #avr
-                      case1234_no_tg, #tg
-                      cases_no_pss) #pss
+gen_OMIB = PSY.DynamicGenerator(
+    1, #Number
+    "OMIB_Gen", #name
+    nodes_OMIB[2], #bus
+    1.0, #ω_ref
+    1.0, #V_ref
+    0.5, #P_ref
+    0.0, #Q_ref: Not used for standard machines (only for PQ gens)
+    machine_OMIB, #machine
+    shaft_OMIB, #shaft
+    avr_OMIB, #avr
+    tg_OMIB, #tg
+    pss_OMIB, #pss
+)
 ```
 
-Note that a generator is defined by its 5 components, while also defining its reference for frequency, voltage and power.
+Note that a generator is defined by its 5 components, while also defining its reference for frequency, voltage and power. The reactive power reference must be defined but is not used for standard machines, since is only used for PQ generators, that will be implemented in future versions of `LITS`.
 
 ### Defining the Dynamic System
 
 Finally, with all the components properly constructed we define the dynamic system:
 
 ```julia
-case1_DynSystem = PSY.System(nodes_case1, #Vector of Buses
-                              branch_case1, #Vector of Branches
-                              [case1_gen], #Vector of Dynamic Injections
-                              vcat(inf_gen_case1,loads_case1), #Vector of Injections
-                              100.0, #MVA Base
-                              60.0) #Freq. Base
+#Create the system
+sys = PSY.System(
+      100.0, #Base MVA
+      frequency = 60.0, #Nominal frequency in Hz
+)
+
+#Add the buses to the system
+for bus in nodes_OMIB
+    PSY.add_component!(sys, bus)
+end
+
+#Add the branches to the system
+for br in branch_OMIB
+    PSY.add_component!(sys, br)
+end
+
+#Add the loads to the system
+for load in loads_OMIB
+    PSY.add_component!(sys, load)
+end
+
+#Add the sources (infinite gens) to the system
+for source in inf_gen_OMIB
+    PSY.add_component!(sys, source)
+end
+
+#Add the generator
+PSY.add_component!(sys, gen_OMIB)
 ```
 
-## Step 3: Initializing the problem
+## Step 3: Build the simulation and initializing the problem
 
-The next step consists in finding an initial condition for the states. But first, we will explore some of the characteristics of our Dynamic System. All information (a ton) can be observed using `dump(case1_DynSystem)`. The following methods can be used to return some information:
-
-- `case1_DynSystem.buses`: Return the vector of buses of the dynamic system.
-- `case1_DynSystem.branches`: Return the vector of branches of the dynamic system.
-- `case1_DynSystem.dyn_injections`: Return the vector of dynamic injections of the dynamic system.
-- `case1_DynSystem.injections`: Return the vector of dynamic injections of the dynamic system.
-- `case1_DynSystem.DAE_vector`: Return the vector of booleans of the dynamic system. Returns `false` or 0 for states that are algebraic and `true` or 1 for states that have derivative defined (differential states). The arrangement will put first the real part of the voltage buses and next the imaginary part. After that the differential states are defined.
-- `case1_DynSystem.global_state_index`: Return an array of dictionaries that have the order of the states in the entire vector state.
-
-To initialize the problem we need to define an initial guess of the states:
+The next step is to create the simulation structure. This will create the indexing of our system that will be used to formulate the differential-algebraic system of equations. To do so, it is required to specify the perturbation that will occur in the system. `LITS` support two types of perturbations:
+- Three Phase Fault
+- Change in Reference Parameter
 
+In here, he will use a Three Phase Fault, that is modeled by modifying the admittance matrix of the system. To do so we create a ThreePhaseFault perturbation as follows:
 ```julia
-#Initialize variables
-dx0 = zeros(LITS.get_total_rows(case1_DynSystem)) #Define a vector of zeros for the derivative
-x0 = [1.05, #VR_1
-      1.0, #VR_2
-      0.0, #VI_1
-      0.01, #VI_2
-      0.2, #δ
-      1.0] #ω
-tspan = (0.0, 30.0);
+#Obtain the Ybus of the faulted system
+Ybus_fault = PSY.Ybus(
+    branch_OMIB_fault, #fault set of lines
+    nodes_OMIB, #set of buses
+)[:,:]
+
+#Construct the perturbation
+perturbation_Ybus = ThreePhaseFault(
+    1.0, #change will occur at t = 1.0s
+    Ybus_fault, #new Ybus
+)
 ```
 
-We will use `NLsolve` to find the initial condition of the system:
-
+With this, we are ready to create our simulation structure. We will skip solving for initial conditions to discuss about indexing. To construct our simulation we use:
 ```julia
-inif! = (out,x) -> LITS.system!(out, #output of the function
-                                      dx0, #derivatives equal to zero
-                                      x, #states
-                                       ([0.0],case1_DynSystem), #Parameters: [0.0] is not used
-                                        0.0) #time equals to zero.
-sys_solve = nlsolve(inif!, x0) #Solve using initial guess x0
-x0_init = sys_solve.zero
+#Time span of our simulation
+tspan = (0.0, 30.0)
+
+#Define Simulation
+sim = Simulation(
+    sys, #system
+    tspan, #time span
+    perturbation_Ybus, #Type of perturbation
+    initialize_simulation = false #keyword argument to not find initial conditions.
+)
 ```
-
-## Step 4: Build the Simulation
-
-Next we will construct the simulation that we are interested to run. But first, we define the pertubation we are interested in model. `LITS` have two perturbations already implemented, that are a change in the mechanical power `P_ref` and a change on the admittance matrix `Y_bus` of the system. In this case we define a change in the admittance matrix:
-
+This will create the simulation structure that will be used to run the transient simulation and will modify the system to include the indexing. `LITS` will have the following structure for the vector of variables:
+```math
+x = \left[\begin{array}{c} v_r \\ v_i \\ z \end{array}\right]
+```
+on which ``v_r`` is the vector of real voltages of all buses, ``v_i`` is the vector of imaginary voltages of all buses and ``z`` is the rest of states defined by the dynamic devices. Then, the length of the vector of variables will ``2n + \text{len}(z)``, where ``n`` is the number of buses in the system. The indexing of the states can be found using:
 ```julia
-#Compute Y_bus after fault
-Ybus_fault = PSY.Ybus(branch_case1_fault, nodes_case1)[:,:] #Obtain Ybus for fault system
-
-#Define Fault using Callbacks
-cb = DiffEqBase.DiscreteCallback(LITS.change_t_one, #Change occurs at t=1
-                                 LITS.Y_change!) #Callback will change the Y_bus.
+ext = PSY.get_ext(sim.system) #Obtain ext information of the system
+ext["global_index"] #Showcase the global indexing of z
 ```
+In this system, ``\delta`` of the generator is state 5 and ``\omega`` is state 6 (since the first 4 states are the bus voltages). In addition, `ext["lits_counts"]` has information on the total variables and total states (differential variables).
 
-Now we define the simulation structure:
+The next step consists in finding an initial condition for the states. In this case simply running
+```julia
+#Define Simulation
+sim = Simulation(
+    sys, #system
+    tspan, #time span
+    perturbation_Ybus, #Type of perturbation
+)
+```
+will correctly initialize the system. If no initial guess is provided, the system will use a flat start guess, assuming that all real voltages are equal to one, while imaginary voltages are equal to zero. Differential variables (states) will be guessed as zero too. The initial values can be obtained using `sim.x0_init`. However, for most systems if a bad initial guess is used, the non-linear solver may fail in correctly initializing the system. For such purposes, an initial guess can be provided to the simulation as follows:
 
 ```julia
-#Define Simulation Problem
-sim = Simulation(case1_DynSystem, #Dynamic System
-                        tspan, #Time span to simulate
-                        Ybus_fault, #Parameter that will be changed in the fault
-                        cb, #Callback
-                        x0_init) #Initial condition
+#Initial guess
+x0_guess = [
+    1.0, #VR_1
+    1.0, #VR_2
+    0.0, #VI_1
+    0.0, #VI_2
+    0.2, #δ
+    1.0, #ω
+]
+
+#Define Simulation
+sim = Simulation(
+    sys, #system
+    tspan, #time span
+    perturbation_Ybus, #Type of perturbation
+    initial_guess = x0_guess, #initial guess
+)
+
+#Check the initial condition
+sim.x0_init
 ```
 
-Finally, to run the simulation:
+## Step 4: Run the Simulation
+
+Finally, to run the simulation we simply use:
 ```julia
-#Solve problem in equilibrium
+#Solve problem
 run_simulation!(sim, #simulation structure
                 IDA(), #Sundials DAE Solver
                 dtmax=0.02); #Arguments: Maximum timestep allowed
 ```
+In some cases, the dynamic time step used for the simulation may fail. In such case, the keyword argument `dtmax` can be used to limit the maximum time step allowed for the simulation.
 
 ## Step 5: Exploring the solution
 
 After running the simulation, our simulation structure `sim` will have the solution. For that `sim.solution` can be used to explore the solution structure. In this case `sim.solution.t` returns the vector of time, while `sim.solution.u` return the array of states. In addition, `LITS` have two functions to obtain different states of the solution:
 
-- `get_state_series(sim, (:Case1Gen, :δ))`: can be used to obtain the solution as a tuple of time and the required state. In this case, we are obtaining the rotor angle `:δ` of the generator named `:Case1Gen`.
+- `get_state_series(sim, ("OMIB_Gen", :δ))`: can be used to obtain the solution as a tuple of time and the required state. In this case, we are obtaining the rotor angle `:δ` of the generator named `"OMIB_Gen"`.
 - `get_voltagemag_series(sim, 2)`: can be used to obtain the voltage magnitude as a tuple of time and voltage. In this case, we are obtaining the voltage magnitude at bus 2 (where the generator is located).
 
 ```julia
 using Plots
-angle = get_state_series(sim, (:Case1Gen, :δ))
+angle = get_state_series(sim, ("OMIB_Gen", :δ))
 plot(angle, xlabel="time", ylabel="rotor angle [rad]", label="rotor angle")
 
 volt = get_voltagemag_series(sim, 2)
@@ -259,3 +320,19 @@ plot(volt, xlabel="time", ylabel="Voltage [pu]", label="V_2")
 ```@raw html
 <img src="../../assets/voltage_OMIB.png" width="75%"/>
 ``` ⠀
+
+## Optional: Small Signal Analysis
+
+`LITS 0.3.0` uses automatic differentiation to compute the reduced Jacobian of the system for the differential states. This can be used to analyze the local stability of the linearized system.
+
+```julia
+small_sig = small_signal_analysis(sim)
+```
+
+The `small_sig` result can report the reduced jacobian for ``\delta`` and ``\omega``, and can also be used to report the eigenvalues of the reduced linearized system.
+
+```julia
+small_sig.reduced_jacobian
+
+small_sig.eigenvalues
+```
diff --git a/docs/src/Examples/example_lines.md b/docs/src/Examples/example_lines.md
index 7e7da1a..c20d688 100644
--- a/docs/src/Examples/example_lines.md
+++ b/docs/src/Examples/example_lines.md
@@ -1,9 +1,11 @@
 # Tutorial: Dynamic Lines
 
 This tutorial will introduce an example of considering dynamic lines in `LITS`.
-Note that this tutorial is for `LITS 0.2.0`. Future versions will have dedicated functions to find an equilibrium point and a proper functions for running perturbations without coding directly the callbacks.
+Note that this tutorial is for `LITS 0.3.0`. Future versions will have dedicated functions to find an equilibrium point using a Power Flow method without relying in a guess of the initial condition to run a non-linear solver.
 
-This tutorial presents a simulation of a three-bus system, with an infinite bus (represented as a voltage source behind an impedance) at bus 1, and a one d- one q- machine on bus 2 and an inverter of 19 states, as a virtual synchronous machine at bus 3. The perturbation will be the trip of two of the three circuits (triplicating its resistance and impedance) of the line that connects bus 1 and bus 3. This case also consider a dynamic line model for connection between buses 2 and 3.
+This tutorial presents a simulation of a three-bus system, with an infinite bus (represented as a voltage source behind an impedance) at bus 1, a one d- one q- machine on bus 2 and an inverter of 19 states, as a virtual synchronous machine at bus 3. The perturbation will be the trip of two of the three circuits (triplicating its resistance and impedance) of the line that connects bus 1 and bus 3. This case also consider a dynamic line model for connection between buses 2 and 3.
+
+It is recommended to check `Tutorial 1: OMIB` first, since that includes more details and explanations on all definitions and functions.
 
 This tutorial can be found on [LITS/Examples](https://github.com/Energy-MAC/LITS-Examples) repository.
 
@@ -12,8 +14,6 @@ This tutorial can be found on [LITS/Examples](https://github.com/Energy-MAC/LITS
 ```julia
 using LITS
 using PowerSystems
-using NLsolve
-using DiffEqBase
 using Sundials
 const PSY = PowerSystems
 ```
@@ -25,34 +25,109 @@ To start we will define the data structures for the network.
 ### Buses and Branches
 
 ```julia
-nodes_case9   = [Bus(1 , "Bus 1"  , "REF" , 0 , 1.02  , (min=0.94, max=1.06), 138),
-                    Bus(2 , "Bus 2"  , "PV" , 0 , 1.00 , (min=0.94, max=1.06), 138),
-                    Bus(3 , "Bus 3"  , "PQ" , 0 , 1.00 , (min=0.94, max=1.06), 138)]
-
-branch_case9  =  [Line("Line1", true, 0.0, 0.0, Arc(from=nodes_case9[1], to=nodes_case9[3]), 0.01, 0.12, (from=0.1, to=0.1), 100, 1.04),
-                    Line("Line2", true, 0.0, 0.0, Arc(from=nodes_case9[1], to=nodes_case9[2]), 0.01, 0.12, (from=0.1, to=0.1), 100, 1.04)]
+nodes_case9 = [
+    PSY.Bus(1, "Bus 1", "REF", 0, 1.02, (min = 0.94, max = 1.06), 138),
+    PSY.Bus(2, "Bus 2", "PV", 0, 1.00, (min = 0.94, max = 1.06), 138),
+    PSY.Bus(3, "Bus 3", "PQ", 0, 1.00, (min = 0.94, max = 1.06), 138),
+]
+
+branch_case9 = [
+    PSY.Line(
+        "Line1",
+        true,
+        0.0,
+        0.0,
+        Arc(from = nodes_case9[1], to = nodes_case9[3]),
+        0.01,
+        0.12,
+        (from = 0.1, to = 0.1),
+        100,
+        1.04,
+    ),
+    PSY.Line(
+        "Line2",
+        true,
+        0.0,
+        0.0,
+        Arc(from = nodes_case9[1], to = nodes_case9[2]),
+        0.01,
+        0.12,
+        (from = 0.1, to = 0.1),
+        100,
+        1.04,
+    ),
+    PSY.Line(
+        "Line3",
+        true,
+        0.0,
+        0.0,
+        Arc(from = nodes_case9[2], to = nodes_case9[3]),
+        0.02,
+        0.9,
+        (from = 0.5, to = 0.5),
+        100,
+        1.04,
+    ),
+]
 
 #Trip of Line 1. Triplicates its impedance
-branch_case9_fault = [Line("Line1", true, 0.0, 0.0, Arc(from=nodes_case9[1], to=nodes_case9[3]), 0.03, 0.36, (from=0.03, to=0.03), 100, 1.04),
-    Line("Line2", true, 0.0, 0.0, Arc(from=nodes_case9[1], to=nodes_case9[2]), 0.01, 0.12, (from=0.1, to=0.1), 100, 1.04)]
-
-#Dynamic Branch between nodes 2 and 3.
-dyn_branch9 = [LITS.DynLine("Line3", true, Arc(from=nodes_case9[2], to=nodes_case9[3]), 0.02, 0.9, (from=0.5, to=0.5))];
+branch_case9_fault = [
+    PSY.Line(
+        "Line1",
+        true,
+        0.0,
+        0.0,
+        Arc(from = nodes_case9[1], to = nodes_case9[3]),
+        0.03,
+        0.36,
+        (from = 0.03, to = 0.03),
+        100,
+        1.04,
+    ),
+    PSY.Line(
+        "Line2",
+        true,
+        0.0,
+        0.0,
+        Arc(from = nodes_case9[1], to = nodes_case9[2]),
+        0.01,
+        0.12,
+        (from = 0.1, to = 0.1),
+        100,
+        1.04,
+    ),
+    PSY.Line(
+        "Line3",
+        true,
+        0.0,
+        0.0,
+        Arc(from = nodes_case9[2], to = nodes_case9[3]),
+        0.02,
+        0.9,
+        (from = 0.5, to = 0.5),
+        100,
+        1.04,
+    ),
+];
 ```
 
 ### Injection devices
 
 ```julia
-loads_case9 =   [PowerLoad("Bus1", true, nodes_case9[1], PowerSystems.ConstantPower, 0.5, 0.1, 1.5, 0.8),
-                   PowerLoad("Bus2", true, nodes_case9[2], PowerSystems.ConstantPower, 1.0, 0.3, 1.5, 0.8),
-                   PowerLoad("Bus3", true, nodes_case9[3], PowerSystems.ConstantPower, 0.3, 0.1, 0.5, 0.3)]
-
-inf_gen_case9 = StaticSource(1, #number
-                    :InfBus, #name
-                    nodes_case9[1],#bus
-                    1.00, #VR
-                    0.0, #VI
-                    0.000005); #Xth
+loads_case9 = [
+    PowerLoad("Bus1", true, nodes_case9[1], PowerSystems.ConstantPower, 0.5, 0.1, 1.5, 0.8),
+    PowerLoad("Bus2", true, nodes_case9[2], PowerSystems.ConstantPower, 1.0, 0.3, 1.5, 0.8),
+    PowerLoad("Bus3", true, nodes_case9[3], PowerSystems.ConstantPower, 0.3, 0.1, 0.5, 0.3),
+]
+
+inf_gen_case9 = [PSY.Source(
+    "InfBus", #name
+    true, #availability
+    nodes_case9[1],#bus
+    1.00, #VR
+    0.0, #VI
+    0.000005,  #Xth
+)]
 ```
 
 ### Dynamic injection devices
@@ -62,201 +137,286 @@ First we define our generator data:
 ```julia
 ######## Machine Data #########
 
-### Case: 4th Order Model with AVR (3-bus case) ###
-case9_machine =  OneDOneQMachine(0.0, #R
-                                  1.3125, #Xd
-                                  1.2578, #Xq
-                                  0.1813, #Xd_p
-                                  0.25, #Xq_p
-                                  5.89, #Td0_p
-                                  0.6, #Tq0_p
-                                  100.0)   #MVABase
+### Case 2: 4th Order Model with AVR (3-bus case) ###
+case9_machine = PSY.OneDOneQMachine(
+    0.0, #R
+    1.3125, #Xd
+    1.2578, #Xq
+    0.1813, #Xd_p
+    0.25, #Xq_p
+    5.89, #Td0_p
+    0.6, #Tq0_p
+    100.0,
+)   #MVABase
 
 ######## Shaft Data #########
 
 ### Shafts for Gen ###
-case9_shaft = SingleMass(3.01, #H (M = 6.02 -> H = M/2)
-                            0.0) #D
+case9_shaft = PSY.SingleMass(
+    3.01, #H (M = 6.02 -> H = M/2)
+    0.0,
+) #D
 
 ######## PSS Data #########
-
-### No PSS ###
-cases_no_pss = PSSFixed(0.0)
-
+cases_no_pss = PSY.PSSFixed(0.0)
 
 ######## TG Data #########
 
-### No TG ###
-case9_no_tg = TGFixed(1.0) #eff
-
+### No TG for Cases 1, 2, 3, 4 ###
+case9_no_tg = PSY.TGFixed(1.0) #eff
 
 ########  AVR Data #########
-### AVRs for this case ###
-case9_avr = AVRTypeI(20.0, #Ka - Gain
-                        0.01, #Ke
-                        0.063, #Kf
-                        0.2, #Ta
-                        0.314, #Te
-                        0.35, #Tf
-                        0.001, #Tr
-                        5.0, #Vrmax
-                        -5.0, #Vrmin
-                        0.0039, #Ae - 1st ceiling coefficient
-                        1.555) #Be - 2nd ceiling coefficient
-
-### Generators ###
-case9_gen = PSY.DynamicGenerator(1, #Number
-                         :Case9Gen,
-                         nodes_case9[2], #bus
-                         1.0, # ω_ref,
-                         1.0124, #V_ref
-                         0.6, #P_ref
-                         case9_machine, #machine
-                         case9_shaft, #shaft
-                         case9_avr, #avr
-                         case9_no_tg, #tg
-                         cases_no_pss); #pss
+### AVRs for Case 2, 3, 4 and 5 ###
+case9_avr = PSY.AVRTypeI(
+    20.0, #Ka - Gain
+    0.01, #Ke
+    0.063, #Kf
+    0.2, #Ta
+    0.314, #Te
+    0.35, #Tf
+    0.001, #Tr
+    5.0, #Vrmax
+    -5.0, #Vrmin
+    0.0039, #Ae - 1st ceiling coefficient
+    1.555,
+) #Be - 2nd ceiling coefficient
+
+### Case 7 Generators ###
+case9_gen = PSY.DynamicGenerator(
+    1, #Number
+    "Case9Gen",
+    nodes_case9[2], #bus
+    1.0, # ω_ref,
+    1.0124, #V_ref
+    0.6, #P_ref
+    0.0, #Q_ref
+    case9_machine, #machine
+    case9_shaft, #shaft
+    case9_avr, #avr
+    case9_no_tg, #tg
+    cases_no_pss, #pss
+);
 ```
 
 and for the inverter:
 
 ```julia
-############### Inverter Data ########################
-
-converter = AvgCnvFixedDC(138.0, #Rated Voltage
-                          100.0) #Rated MVA
-
-dc_source = FixedDCSource(1500.0) #Not in the original data, guessed.
-
-filt = LCLFilter(0.08, #Series inductance lf in pu
-                   0.003, #Series resitance rf in pu
-                   0.074, #Shunt capacitance cf in pu
-                   0.2, #Series reactance rg to grid connection (#Step up transformer or similar)
-                   0.01) #Series resistance lg to grid connection (#Step up transformer or similar)
-
-pll = PLL(500.0, #ω_lp: Cut-off frequency for LowPass filter of PLL filter.
-          0.084, #k_p: PLL proportional gain
-          4.69) #k_i: PLL integral gain
-
-virtual_H = VirtualInertia(2.0, #Ta:: VSM inertia constant
-                           400.0, #kd:: VSM damping coefficient
-                           20.0, #kω:: Frequency droop gain in pu
-                           2*pi*50.0) #ωb:: Rated angular frequency
-
-Q_control = ReactivePowerDroop(0.2, #kq:: Reactive power droop gain in pu
-                              1000.0) #ωf:: Reactive power cut-off low pass filter frequency
-
-outer_control = VirtualInertiaQdroop(virtual_H, Q_control)
-
-vsc = CombinedVIwithVZ(0.59, #kpv:: Voltage controller proportional gain
-                       736.0, #kiv:: Voltage controller integral gain
-                       0.0, #kffv:: Binary variable enabling the voltage feed-forward in output of current controllers
-                       0.0, #rv:: Virtual resistance in pu
-                       0.2, #lv: Virtual inductance in pu
-                       1.27, #kpc:: Current controller proportional gain
-                       14.3, #kiv:: Current controller integral gain
-                       0.0, #kffi:: Binary variable enabling the current feed-forward in output of current controllers
-                       50.0, #ωad:: Active damping low pass filter cut-off frequency
-                       0.2) #kad:: Active damping gain
-
-case9_inv = PSY.DynamicInverter(2, #number
-                             :DARCO, #name
-                             nodes_case9[3], #bus location
-                             1.0, #ω_ref
-                             0.8, #V_ref
-                             0.5, #P_ref
-                             -0.3, #Q_ref
-                             100.0, #MVABase
-                             converter, #Converter
-                             outer_control, #OuterControl
-                             vsc, #Voltage Source Controller
-                             dc_source, #DC Source
-                             pll, #Frequency Estimator
-                             filt); #Output Filter
+converter = PSY.AvgCnvFixedDC(
+    138.0, #Rated Voltage
+    100.0,
+) #Rated MVA
+
+dc_source = PSY.FixedDCSource(1500.0) #Not in the original data, guessed.
+
+filt = PSY.LCLFilter(
+    0.08, #Series inductance lf in pu
+    0.003, #Series resitance rf in pu
+    0.074, #Shunt capacitance cf in pu
+    0.2, #Series reactance rg to grid connection (#Step up transformer or similar)
+    0.01,
+) #Series resistance lg to grid connection (#Step up transformer or similar)
+
+pll = PSY.PLL(
+    500.0, #ω_lp: Cut-off frequency for LowPass filter of PLL filter.
+    0.084, #k_p: PLL proportional gain
+    4.69,
+) #k_i: PLL integral gain
+
+virtual_H = PSY.VirtualInertia(
+    2.0, #Ta:: VSM inertia constant
+    400.0, #kd:: VSM damping coefficient
+    20.0, #kω:: Frequency droop gain in pu
+    2 * pi * 50.0,
+) #ωb:: Rated angular frequency
+
+Q_control = PSY.ReactivePowerDroop(
+    0.2, #kq:: Reactive power droop gain in pu
+    1000.0,
+) #ωf:: Reactive power cut-off low pass filter frequency
+
+outer_control = PSY.VirtualInertiaQdroop(virtual_H, Q_control)
+
+vsc = PSY.CombinedVIwithVZ(
+    0.59, #kpv:: Voltage controller proportional gain
+    736.0, #kiv:: Voltage controller integral gain
+    0.0, #kffv:: Binary variable enabling the voltage feed-forward in output of current controllers
+    0.0, #rv:: Virtual resistance in pu
+    0.2, #lv: Virtual inductance in pu
+    1.27, #kpc:: Current controller proportional gain
+    14.3, #kiv:: Current controller integral gain
+    0.0, #kffi:: Binary variable enabling the current feed-forward in output of current controllers
+    50.0, #ωad:: Active damping low pass filter cut-off frequency
+    0.2,
+) #kad:: Active damping gain
+
+case9_inv = PSY.DynamicInverter(
+    2, #number
+    "DARCO", #name
+    nodes_case9[3], #bus location
+    1.0, #ω_ref
+    0.8, #V_ref
+    0.5, #P_ref
+    -0.3, #Q_ref
+    100.0, #MVABase
+    converter, #Converter
+    outer_control, #OuterControl
+    vsc, #Voltage Source Controller
+    dc_source, #DC Source
+    pll, #Frequency Estimator
+    filt, #Output Filter
+);
 ```
 
-### Dynamic System
+### Defining the Dynamic System
 
 ```julia
-case9_DynSystem = PSY.System(nodes_case9, #Vector of Buses
-                                 branch_case9, #Vector of Branches
-                                 [case9_inv, case9_gen], #Vector of Dynamic Injections
-                                 vcat(inf_gen_case9, loads_case9), #Vector of Injections
-                                 100.0, #MVA Base
-                                 50.0, #Freq. Base
-                                 dyn_branch9); #Vector of Dynamic Branches
+#Create the system
+sys = PSY.System(
+      100.0, #Base MVA
+      frequency = 50.0, #Nominal frequency in Hz
+)
+
+#Add the buses to the system
+for bus in nodes_case9
+    PSY.add_component!(sys, bus)
+end
+
+#Add the branches to the system
+for br in branch_case9
+    PSY.add_component!(sys, br)
+end
+
+#Make Line3 a dynamic line
+make_dynamic_branch!(branch_case9[3], sys)
+
+#Add the loads to the system
+for load in loads_case9
+    PSY.add_component!(sys, load)
+end
+
+#Add the sources (infinite gens) to the system
+for source in inf_gen_case9
+    PSY.add_component!(sys, source)
+end
+
+#Add the inverter
+PSY.add_component!(sys, case9_inv)
+
+#Add the generator
+PSY.add_component!(sys, case9_gen)
 ```
 
-## Step 3: Initializing the problem
+## Step 3: Build the simulation and initializing the problem
 
-To initialize the problem we need to define an initial guess of the states:
+First, we construct the perturbation, by properly computing the new Ybus:
 
 ```julia
-dx0 = zeros(LITS.get_total_rows(case9_DynSystem))
-x0 = [1.00, #V1_R
-          1.00, #V2_R
-          1.00, #V3_R
-          0.0, #V1_I
-         -0.01, #V2_I
-         -0.01, #V3_I
-         0.0, #δω_vsm
-          0.2, #δθ_vsm
-          0.025, #qm
-          0.0015, #ξ_d
-          -0.07, #ξ_q
-          0.05, #γ_d
-        -0.001, #γ_q
-         0.95, #ϕ_d
-         -0.10, #ϕ_q
-         1.004, #vpll_d
-         0.0, #vpll_q
-         0.0, #ε_pll
-         0.1, #δθ_pll
-         0.5, #id_cv
-         0.0, #iq_cv
-         0.95, #vod
-         -0.1, #voq
-         0.49, #iod
-        -0.1, #ioq
-          1.0, #eq_p
-          0.47, #ed_p
-          0.6, #δ
-          1.0, #ω
-          2.1, #Vf
-          0.28, #Vr1
-          -0.39, #Vr2,
-          1.0, #Vm
-          0.5, #IL1_R
-          0.5] #IL1_I
-
-case9_inv.inner_vars[13] = 0.95 #Initialize internal variables of inverter: Vd_cnv var
-case9_inv.inner_vars[14] = -0.1 #Initialize internal variables of inverter: Vq_cnv var
-tspan = (0.0, 30.0);
+#Obtain perturbed Ybus
+sys2 = PSY.System(
+      100.0, #Base MVA
+      frequency = 50.0, #Nominal frequency in Hz
+)
+
+#Add the buses to the system
+for bus in nodes_case9
+    PSY.add_component!(sys2, bus)
+end
+
+#Add the fault branches to the system
+for br in branch_case9_fault
+    PSY.add_component!(sys2, br)
+end
+
+#Make Line3 a dynamic line
+make_dynamic_branch!(branch_case9_fault[3], sys2)
+
+#Compute the Ybus
+Ybus_fault = PSY.Ybus(sys2)[:,:]
+
+#Construct the perturbation
+perturbation_Ybus = ThreePhaseFault(
+    1.0, #change will occur at t = 1.0s
+    Ybus_fault, #new Ybus
+);
 ```
 
-We will use `NLsolve` to find the initial condition of the system:
+Now, we check the indexing
 
 ```julia
-#Find initial condition
-inif! = (out,x) -> LITS.system!(out, dx0 ,x, ([0.0],case9_DynSystem), 0.0)
-sys_solve = nlsolve(inif!, x0, xtol=:1e-8,ftol=:1e-8,method=:trust_region)
-x0_init = sys_solve.zero
+#Time span of our simulation
+tspan = (0.0, 30.0)
+
+#Define Simulation
+sim = Simulation(
+    sys, #system
+    tspan, #time span
+    perturbation_Ybus, #Type of perturbation
+    initialize_simulation = false #keyword argument to not find initial conditions.
+)
+
+ext = PSY.get_ext(sim.system) #Obtain ext information of the system
+ext["global_index"] #Showcase the global indexing of z
 ```
 
-## Step 4: Build the Simulation
+This is a clear example on which a flat start is not enough to properly initialize the system, so we provide a guess:
 
 ```julia
-#Compute Y_bus after fault
-Ybus_fault = Ybus(branch_case9_fault, nodes_case9)[:,:]
+x0_guess = [
+    1.00, #V1_R
+    1.00, #V2_R
+    1.00, #V3_R
+    0.0, #V1_I
+    -0.01, #V2_I
+    -0.01, #V3_I
+    0.0, #δω_vsm
+    0.2, #δθ_vsm
+    0.025, #qm
+    0.0015, #ξ_d
+    -0.07, #ξ_q
+    0.05, #γ_d
+    -0.001, #γ_q
+    0.95, #ϕ_d
+    -0.10, #ϕ_q
+    1.004, #vpll_d
+    0.0, #vpll_q
+    0.0, #ε_pll
+    0.1, #δθ_pll
+    0.5, #id_cv
+    0.0, #iq_cv
+    0.95, #vod
+    -0.1, #voq
+    0.49, #iod
+    -0.1, #ioq
+    1.0, #eq_p
+    0.47, #ed_p
+    0.6, #δ
+    1.0, #ω
+    2.1, #Vf
+    0.28, #Vr1
+    -0.39, #Vr2,
+    1.0, #Vm
+    0.5, #IL1_R
+    0.5, #IL1_I
+]
+```
+and with that we can compute our initial condition:
+```julia
+sim = Simulation(
+    sys, #system
+    tspan, #time span
+    perturbation_Ybus, #Type of perturbation
+    initial_guess = x0_guess #initial guess.
+)
+```
 
-#Define Fault using Callbacks
-cb = DiffEqBase.DiscreteCallback(LITS.change_t_one, LITS.Y_change!)
 
-#Define Simulation Problem
-sim = Simulation(case9_DynSystem, tspan, Ybus_fault, cb, x0_init)
+## Step 4: Run the simulation
 
-#Solve problem
-run_simulation!(sim, IDA());
+```julia
+#Run the simulation
+run_simulation!(sim, #simulation structure
+                IDA(), #Sundials DAE Solver
+)
 ```
 
 ## Step 5: Explore the solution
diff --git a/docs/src/Models/gens.md b/docs/src/Models/gens.md
index feaa053..13164cd 100644
--- a/docs/src/Models/gens.md
+++ b/docs/src/Models/gens.md
@@ -22,13 +22,13 @@ Each generator is defined in its own ``dq`` reference frame. Let ``\delta`` be t
 \end{align}
 ```
 
-Models are based from Federico Milano's book: "Power System Modelling and Scripting" and Prabha Kundur's book: "Power System's Stability and Control"
+Models are based from Federico Milano's book: "Power System Modelling and Scripting" and Prabha Kundur's book: "Power System's Stability and Control" and structures are defined in ```PowerSystems.jl``` abbreviated as ```PSY```.
 
 ## Machines
 
 The machine component describes the stator-rotor electromagnetic dynamics.
 
-### Classical Model (Zero Order)
+### Classical Model (Zero Order) ```[PSY.BaseMachine]```
 This is the classical order model that does not have differential equations in its machine model (``\delta`` and ``\omega`` are defined in the shaft):
 
 ```math
@@ -38,7 +38,7 @@ p_e \approx \tau_e &= (v_q + r_a i_q)i_q + (v_d + r_ai_d)i_d \tag{1b}
 \end{align}
 ```
 
-### One d- One q- Model (2th Order)
+### One d- One q- Model (2th Order) ```[PSY.OneDOneQMachine]```
 This model includes two transient emf with their respective differential equations:
 
 ```math
@@ -50,7 +50,7 @@ p_e \approx \tau_e &= (v_q + r_a i_q)i_q + (v_d + r_ai_d)i_d \tag{2d}
 \end{align}
 ```
 
-### Marconato Machine (6th Order)
+### Marconato Machine (6th Order) ```[PSY.MarconatoMachine]```
 
 The Marconato model defines 6 differential equations, two for stator fluxes and 4 for transient and subtransient emfs:
 
@@ -77,7 +77,7 @@ with
 \end{align*}
 ```
 
-### Simplified Marconato Machine (4th Order)
+### Simplified Marconato Machine (4th Order) ```[PSY.SimpleMarconatoMachine]```
 
 This model neglects the derivative of stator fluxes (``\dot{\psi}_d`` and  ``\dot{\psi}_q``) and assume that the rotor speed stays close to 1 pu (``\omega\psi_{d}=\psi_{d}`` and ``\omega\psi_{q}=\psi_{q}``) that allows to remove the stator fluxes variables from the Marconato model.
 
@@ -102,7 +102,7 @@ with
 ```
 
 
-### Anderson-Fouad Machine (6th Order)
+### Anderson-Fouad Machine (6th Order) ```[PSY.AndersonFouadMachine]```
 
 The Anderson-Fouad model also defines 6 differential equations, two for stator fluxes and 4 for transient and subtransient emfs and is derived from the Marconato model by defining ``\gamma_d \approx \gamma_q \approx T_{AA} \approx 0``:
 
@@ -120,7 +120,7 @@ i_q &= \frac{1}{x_q''} (-e_d'' - \psi_q) \tag{5h} \\
 \end{align}
 ```
 
-### Simplified Anderson-Fouad Machine (4th Order)
+### Simplified Anderson-Fouad Machine (4th Order) ``[PSY.SimpleAFMachine]``
 
 Similar to the Simplified Marconato Model, this model neglects the derivative of stator fluxes (``\dot{\psi}_d`` and  ``\dot{\psi}_q``) and assume that the rotor speed stays close to 1 pu (``\omega \psi_d = \psi_d`` and ``\omega \psi_q = \psi_q``) that allows to remove the stator fluxes variables from the model:
 
@@ -135,7 +135,7 @@ p_e \approx \tau_e &= (v_q + r_a i_q)i_q + (v_d + r_ai_d)i_d \tag{6f}
 \end{align}
 ```
 
-### Rotor Fluxes Model (5th Order)
+### Rotor Fluxes Model (5th Order) ```[PSY.FullMachine]```
 
 This models describes a machine model using 2 stator fluxes ``\psi_d`` and ``\psi_q``, the rotor field flux ``\psi_f``, the d-axis damping ``\psi_{1d}`` and only one q-axis damping circuit ``\psi_{1q}``:
 
@@ -157,7 +157,7 @@ This models describes a machine model using 2 stator fluxes ``\psi_d`` and ``\ps
 
 The shaft component defines the rotating mass of the synchronous generator.
 
-### Rotor Mass Shaft
+### Rotor Mass Shaft ```[PSY.SingleMass]```
 
 This is the standard model, on which one single mass (typically the rotor) is used to model the entire inertia of the synchronous generator. Each generator's rotating frame use a reference frequency ``\omega_s``, that typically is the synchronous one (i.e. ``\omega_s = 1.0``). The model defines two differential equations for the rotor angle ``\delta`` and the rotor speed ``\omega``:
 
@@ -169,7 +169,7 @@ This is the standard model, on which one single mass (typically the rotor) is us
 ```
 
 
-### Five-Mass Shaft
+### Five-Mass Shaft ```[PSY.FiveMassShaft]```
 
 This model describes model connecting a high-pressure (hp) steam turbine, intermediate-pressure (ip) steam turbine, low-pressure (lp) steam pressure, rotor and exciter (ex) connected in series (in that order) in the same shaft using a spring-mass model:
 
@@ -192,11 +192,11 @@ This model describes model connecting a high-pressure (hp) steam turbine, interm
 
 AVR are used to determine the voltage in the field winding ``v_f`` in the model.
 
-### Fixed AVR
+### Fixed AVR ```[PSY.AVRFixed]```
 
 This is a simple model that set the field voltage to be equal to a desired constant value ``v_f = v_{\text{fix}}``.
 
-### Simple AVR
+### Simple AVR ```[PSY.AVRSimple]```
 
 This depicts the most basic AVR, on which the field voltage is an integrator over the difference of the measured voltage and a reference:
 
@@ -206,7 +206,7 @@ This depicts the most basic AVR, on which the field voltage is an integrator ove
 \end{align}
 ```
 
-### AVR Type I
+### AVR Type I ```[PSY.AVRTypeI]```
 
 This AVR is a simplified version of the IEEE DC1 AVR model:
 
@@ -227,7 +227,7 @@ S_e(v_f) = A_e \exp(B_e|v_f|)
 \end{align*}
 ```
 
-### AVR Type II
+### AVR Type II ```[PSY.AVRTypeII]```
 
 This model represents a static exciter with higher gains and faster response than the Type I:
 
@@ -254,11 +254,11 @@ S_e(v_f) &= A_e \exp(B_e|v_f|)
 
 PSS are used to add an additional signal ``v_s`` to the field voltage: ``v_f = v_f^{\text{avr}} + v_s``.
 
-### Fixed PSS
+### Fixed PSS ```[PSY.PSSFixed]```
 
 This is a simple model that set the stabilization signal to be equal to a desired constant value ``v_s = v_{s}^{\text{fix}}``. The absence of PSS can be modelled using this component with ``v_s^{\text{fix}} = 0``.
 
-### Simple PSS
+### Simple PSS ```[PSY.PSSSimple]```
 
 This is the most basic PSS that can be implemented, on which the stabilization signal is simply a proportional controller over the frequency and electrical power:
 
@@ -272,11 +272,11 @@ v_s = K_{\omega}(\omega - \omega_s) + K_p(\omega \tau_e - P_{\text{ref}}) \tag{1
 
 This section describes how mechanical power is modified to provide primary frequency control with synchronous generators. It is assumed that ``\tau_{\text{ref}} = P_{\text{ref}}`` since they are decided at nominal frequency ``\omega = 1``.
 
-### Fixed TG
+### Fixed TG ```[PSY.TGFixed]```
 
 This a simple model that set the mechanical torque to be equal to a proportion of the desired reference ``\tau_m = \eta P_{\text{ref}}``. To set the mechanical torque to be equal to the desired power, the value of ``\eta`` is set to 1.
 
-### TG Type I
+### TG Type I ```[PSY.TGTypeI]```
 
 This turbine governor is described by a droop controller and a low-pass filter to model the governor and two lead-lag blocks to model the servo and reheat of the turbine governor.
 
@@ -297,7 +297,7 @@ p_{\text{in}} = P_{\text{ref}} + \frac{1}{R}(\omega_s - 1)
 \end{align*}
 ```
 
-### TG Type II
+### TG Type II ```[PSY.TGTypeII]```
 
 This turbine governor is a simplified model of the Type I.
 
@@ -311,23 +311,4 @@ This turbine governor is a simplified model of the Type I.
 
 ## Reference
 
-```@docs
-LITS.BaseMachine
-LITS.OneDOneQMachine
-LITS.MarconatoMachine
-LITS.SimpleMarconatoMachine
-LITS.AndersonFouadMachine
-LITS.SimpleAFMachine
-LITS.FullMachine
-LITS.SingleMass
-LITS.FiveMassShaft
-LITS.AVRFixed
-LITS.AVRSimple
-LITS.AVRTypeI
-LITS.AVRTypeII
-LITS.PSSFixed
-LITS.PSSSimple
-LITS.TGFixed
-LITS.TGTypeI
-LITS.TGTypeII
-```
+For constructors check the API on [PowerSystems.jl documentation](https://nrel.github.io/PowerSystems.jl/latest/api/PowerSystems/)
diff --git a/docs/src/Models/inverters.md b/docs/src/Models/inverters.md
index 80b4740..eeed09f 100644
--- a/docs/src/Models/inverters.md
+++ b/docs/src/Models/inverters.md
@@ -17,13 +17,13 @@ The following figure summarizes the components of a inverter and which variables
 
 Contrary to the generator, there are many control structures that can be used to model inverter controllers (e.g. grid-following, grid feeding or virtual synchronous machine). For this purpose, more variables are shared among the components in order to cover all these posibilities.
 
-Models are based from the paper: "A Virtual Synchronous Machine implementation for distributed control of power converters in SmartGrids" from S. D'Arco, J.A. Suul and O.B. Fosso.
+Models are based from the paper: "A Virtual Synchronous Machine implementation for distributed control of power converters in SmartGrids" from S. D'Arco, J.A. Suul and O.B. Fosso, and structures are defined in ```PowerSystems.jl``` abbreviated as ```PSY```.
 
 ## DC Source
 
 This component can be used to model the dynamics of the DC side of the converter.
 
-### Fixed DC Source
+### Fixed DC Source ```[PSY.FixedDCSource]```
 
 This is a model that set the DC voltage to a fixed value ``v_{\text{dc}} = v_{\text{dc}}^{\text{fix}}``.
 
@@ -32,7 +32,7 @@ This is a model that set the DC voltage to a fixed value ``v_{\text{dc}} = v_{\t
 
 This component is used to estimate the frequency of the grid based on the voltage at the bus.
 
-### Phase-Locked Loop (PLL) for VSM
+### Phase-Locked Loop (PLL) for VSM ```[PSY.PLL]```
 
 The following equations present a PLL used to estimate the frequency and PLL angle of the grid. There are two reference frames considered in this inverter. Those are the VSM of the outer-loop control ``\delta\theta_{\text{olc}}`` and the PLL one ``\delta\theta_{\text{pll}}``. The notation used a ``\delta\theta`` to refer as the variation of the respective angle ``\theta`` with respect to the grid SRF (instead of the fixed ``\alpha`` component of the ``\alpha\beta`` transformation):
 
@@ -57,13 +57,13 @@ with
 
 This component defines controllers for both active and reactive power
 
-### Virtual Inertia and Q-droop
+### Virtual Inertia and Q-droop ```[PSY.OuterControl]```
 
-The following model represent a virtual synchronous machine model to represent how active power is going to be deployed. It defines a new SRF denoted as ``\theta_{\text{olc}}`` for the active power controller and uses a simple voltage droop for dispatching reactive power:
+The following model represent a virtual synchronous machine model to represent how active power is going to be deployed. The constructor is ```PSY.OuterControl{PSY.VirtualInertia, PSY.ReactivePowerDroop}```. It defines a new SRF denoted as ``\theta_{\text{olc}}`` for the active power controller and uses a simple voltage droop for dispatching reactive power:
 
 ```math
 \begin{align}
-    \dot{\delta\omega}_{\text{olc}} &= \frac{p_{\text{ref}}}{T_a} - \frac{p_e}{T_a} - \frac{k_d(\omega_{\text{vsm}} - \omega_{\text{pll}})}{T_a} - \frac{k_\omega(\omega_{\text{olc}} - \omega_{\text{ref}})}{T_a} \tag{2a} \\
+    \dot{\delta\omega}_{\text{olc}} &= \frac{p_{\text{ref}}}{T_a} - \frac{p_e}{T_a} - \frac{k_d(\omega_{\text{olc}} - \omega_{\text{pll}})}{T_a} - \frac{k_\omega(\omega_{\text{olc}} - \omega_{\text{ref}})}{T_a} \tag{2a} \\
     \dot{\delta\theta}_{\text{olc}} &= \Omega_b \delta\omega_{\text{olc}} \tag{2b} \\
     \dot{q}_m &= \omega_f (q_e - q_m) \tag{2c}
 \end{align}
@@ -83,7 +83,7 @@ with
 
 This component defines voltage and current controllers to generate the reference signal for the converter.
 
-### Integrated Virtual Impedance, Voltage and Current Controller
+### Integrated Virtual Impedance, Voltage and Current Controller ```[PSY.CombinedVIwithVZ]```
 
 The following model receives both the outer loop control frequency and reference voltage signal to generate the reference signal for the converters. The virtual impedance plays a similar role of the impedance of a synchronous generator. A PI voltage controller is used to generate the current signal that is used in the PI current controller to finally generate the voltage reference signal for the converters.
 
@@ -115,7 +115,7 @@ with
 
 This component can be used to model the dynamics of the switching process.
 
-### Average Model
+### Average Model ```[PSY.AvgCnvFixedDC]```
 
 The average model simply output the desired reference signal since:
 
@@ -130,7 +130,7 @@ where ``m_{dq}`` is the modulation signal, and ``v_{dq}^{\text{ref-signal}}`` is
 
 ## Filters
 
-### LCL Filter
+### LCL Filter ```[PSY.LCLFilter]```
 
 A standard LCL filter is proposed to connect the output of the converter to the grid. In this case, ``v_d`` and ``v_q`` are voltages in the capacitor, while ``v_d^{\text{grid}}`` and ``v_q^{\text{grid}}`` represent the voltage at the bus. The L filter after the capacitor can also include a step-up transformer to increase the voltage, that is model as an extra impedance.
 
@@ -148,13 +148,4 @@ A standard LCL filter is proposed to connect the output of the converter to the
 
 ## Reference
 
-```@docs
-LITS.FixedDCSource
-LITS.PLL
-LITS.VirtualInertia
-LITS.ReactivePowerDroop
-LITS.VirtualInertiaQdroop
-LITS.CombinedVIwithVZ
-LITS.AvgCnvFixedDC
-LITS.LCLFilter
-```
+For constructors check the API on [PowerSystems.jl documentation](https://nrel.github.io/PowerSystems.jl/latest/api/PowerSystems/)
diff --git a/docs/src/Models/network.md b/docs/src/Models/network.md
index c547876..a8f857a 100644
--- a/docs/src/Models/network.md
+++ b/docs/src/Models/network.md
@@ -22,10 +22,10 @@ Each line is defined using a ``\pi`` model connecting two buses ``(n,m)``, with
 
 ```math
 \begin{align}
-Y_{nn} &+= \frac{1}{r+jx} + jc_n \\
-Y_{nm} &+= \frac{-1}{r+jx} \\
-Y_{mm} &+= \frac{1}{r+jx} + jc_m \\
-Y_{mn} &+= \frac{-1}{r+jx} \\
+Y_{nn} &+\!= \frac{1}{r+jx} + jc_n \\
+Y_{nm} &+\!= \frac{-1}{r+jx} \\
+Y_{mm} &+\!= \frac{1}{r+jx} + jc_m \\
+Y_{mn} &+\!= \frac{-1}{r+jx} \\
 \end{align}
 ```
 
diff --git a/docs/src/Models/small.md b/docs/src/Models/small.md
new file mode 100644
index 0000000..0c57ab5
--- /dev/null
+++ b/docs/src/Models/small.md
@@ -0,0 +1,62 @@
+# Small Signal Analysis
+
+Here we discuss the method used to do a small signal analysis on the DAE system defined in *LITS.jl*. The package defines algebraic variables for both real and imaginary voltages on all buses (except if they have a dynamic line connected, on which the voltage of those buses are treated as differential variables). In addition, each dynamic device can add differential variables (or states) that are concatenated to construct the system of differential algebraic equations.
+
+We define ``y`` as the vector of algebraic variables, ``x`` as the vector of differential variables (states) and ``p`` the parameters of the system, we can define ``g(y,x,p)`` as the vector of algebraic equations and ``f(y,x,p)`` as the vector of differential equations. With that, the non-linear differential algebraic system of equations can be written as:
+
+```math
+\begin{align}
+\left[\begin{array}{c}
+ 0 \\
+  \dot{x}
+  \end{array}\right] = \left[\begin{array}{c}
+  g(y,x,p) \\
+   f(y,x,p) \end{array}\right]
+\end{align}
+```
+
+For small signal analysis, we are interested in the stability around an equilbrium point ``y_{eq},x_{eq}`` that satisfies ``\dot{x} = 0`` or equivalently ``f(y_{eq},x_{eq},p) = 0``, while obviously satisfying ``g(y_{eq}, x_{eq}, p) = 0``. To do that we use a first order approximation:
+
+```math
+\begin{align}
+\left[\begin{array}{c}
+ 0 \\
+  \Delta\dot{x}
+  \end{array}\right] = \underbrace{\left[\begin{array}
+  ~g(y_{eq},x_{eq},p) \\
+   f(y_{eq},x_{eq},p) \end{array}\right]}_{ =~ 0}
+ + J[y_{eq}, x_{eq}, p] \left[\begin{array}{c}
+ \Delta y \\
+  \Delta x
+  \end{array}\right]
+  \end{align}
+```
+
+The first to note is that the jacobian matrix can be splitted in 4 blocks depending on the specific variables we are taking the partial derivatives:
+```math
+\begin{align}
+J[y_{eq}, x_{eq}, p] =
+\left[\begin{array}{cc}
+ g_y & g_x \\
+ f_y & f_x \\
+  \end{array}\right]
+\end{align}
+```
+For small signal analyses, we are interested in the stability of the differential states, while still considering that those need to evolve in the manifold defined by the linearized algebraic equations. Assuming that ``g_y`` is not singular (see chapter 7 of Federico Milano's book: "Power System Modelling and Scripting" or [the following paper](https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=1323205)) we can eliminate the algebraic variables to obtain the reduced jacobian:
+
+```math
+\begin{align}
+J_{\text{red}} = f_x - f_y g_y^{-1} g_x
+\end{align}
+```
+that defines our reduced system for the differential variables
+```math
+\begin{align}
+\Delta \dot{x} = J_{\text{red}} \Delta x
+\end{align}
+```
+on which we can compute its eigenvalues to analyze local stability.
+
+## Automatic Differentiation
+
+Once an equilibrium point is found, the complete jacobian of the non-linear system can be obtained using [automatic differentiation](https://en.wikipedia.org/wiki/Automatic_differentiation) in [Julia](https://www.juliadiff.org). In particular, the package `ForwardDiff.jl` is used to obtain the jacobian of the non-linear algebraic system of equations. *LITS.jl* handles the resulting jacobian and reports the reduced jacobian and the corresponding eigenvalues and eigenvectors.
diff --git a/docs/src/assets/SoftwareLoop.png b/docs/src/assets/SoftwareLoop.png
index b4c3f7e..b1af4fd 100644
Binary files a/docs/src/assets/SoftwareLoop.png and b/docs/src/assets/SoftwareLoop.png differ
diff --git a/docs/src/index.md b/docs/src/index.md
index ec72f2d..9340b05 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -22,10 +22,10 @@ using LITS
 ## Structure
 
 The following figure shows the interactions between `LITS.jl`, `PowerSystems.jl`, `DifferentialEquations.jl` and the integrators.
-The architecture of `LITS.jl`  is such that the power system models are all self-contained and return the model function evaluations. The Jacobian is calculated through `DifferentialEquations.jl`'s common-interface enabling the use of any solver available in Julia. Considering that the resulting models are differential-algebraic equations (DAE), the implementation focuses on the use of implicit solvers, in particular SUNDIALS since it has exceptional features applicable to large models — for instance, interfacing with distributed linear-solvers and GPU arrays.
+The architecture of `LITS.jl`  is such that the power system models are all self-contained and return the model function evaluations. The Jacobian is calculated through `DifferentialEquations.jl`'s common-interface enabling the use of any solver available in Julia. Considering that the resulting models are differential-algebraic equations (DAE), the implementation focuses on the use of implicit solvers, in particular SUNDIALS since it has exceptional features applicable to large models — for instance, interfacing with distributed linear-solvers and GPU arrays. In addition, automatic differentiation is implemented using `ForwardDiff.jl` to obtain jacobians to perform small signal analysis.
 
 ```@raw html
-<img src="./assets/SoftwareLoop.png" width="50%"/>
+<img src="./assets/SoftwareLoop.png" width="60%"/>
 ``` ⠀
 
 ## Contents
@@ -37,5 +37,6 @@ Pages = [
   "Models/network.md",
   "Models/gens.md",
   "Models/inverters.md",
+  "Models/small.md",
 ]
 ```
diff --git a/src/base/simulation.jl b/src/base/simulation.jl
index 2e50014..2ec9b61 100644
--- a/src/base/simulation.jl
+++ b/src/base/simulation.jl
@@ -25,7 +25,10 @@ function Simulation(
 
     if initialize_simulation
         @info("Initializing Simulation States")
-        x0_guess = get(kwargs, :initial_guess, zeros(var_count))
+        flat_start = zeros(var_count)
+        bus_count = length(PSY.get_components(PSY.Bus, system))
+        flat_start[1:bus_count] .= 1.0
+        x0_guess = get(kwargs, :initial_guess, flat_start)
         x0_init, initialized = _calculate_initial_conditions(system, x0_guess)
     else
         x0_init = zeros(var_count)
diff --git a/test/data_tests/dynamic_test_data.jl b/test/data_tests/dynamic_test_data.jl
index 6f255b6..12f6e4a 100644
--- a/test/data_tests/dynamic_test_data.jl
+++ b/test/data_tests/dynamic_test_data.jl
@@ -373,6 +373,23 @@ function dyn_gen_case8(nodes)
     ) #pss
 end
 
+function dyn_gen_case9(nodes)
+    return PSY.DynamicGenerator(
+        1, #Number
+        "Case9Gen",
+        nodes[2], #bus
+        1.0, # ω_ref,
+        1.0124, #V_ref
+        0.6, #P_ref
+        0.0, #Q_ref
+        machine_4th(), #machine
+        shaft_no_damping(), #shaft
+        avr_type1(), #avr
+        tg_none(), #tg
+        pss_none(),
+    ) #pss
+end
+
 ######################################
 ############# Inverters ##############
 ######################################
@@ -582,3 +599,43 @@ function system_50Hz(nodes, branches, loads, sources, invs, gens)
 
     return sys
 end
+
+function system_case9(nodes, branches, loads, sources, invs, gens)
+    #Create system with BasePower = 100 MVA and nominal frequency 50 Hz.
+    sys = PSY.System(100.0, frequency = 50.0)
+
+    #Add buses
+    for bus in nodes
+        PSY.add_component!(sys, bus)
+    end
+
+    #Add lines
+    for lines in branches
+        PSY.add_component!(sys, lines)
+    end
+
+    #Make line 3 the dynamic line
+    make_dynamic_branch!(branches[3], sys)
+
+    #Add loads
+    for load in loads
+        PSY.add_component!(sys, load)
+    end
+
+    #Add infinite source
+    for source in sources
+        PSY.add_component!(sys, source)
+    end
+
+    #Add inverters
+    for inv in invs
+        PSY.add_component!(sys, inv)
+    end
+
+    #Add generators
+    for gen in gens
+        PSY.add_component!(sys, gen)
+    end
+
+    return sys
+end
diff --git a/test/data_tests/network_test_data.jl b/test/data_tests/network_test_data.jl
index 9b3acb3..9ce3267 100644
--- a/test/data_tests/network_test_data.jl
+++ b/test/data_tests/network_test_data.jl
@@ -304,6 +304,84 @@ branches_3lines_case8_fault(nodes_3bus) = [
     ),
 ]
 
+branch_3bus_case9(nodes) = [
+    PSY.Line(
+        "Line1",
+        true,
+        0.0,
+        0.0,
+        Arc(from = nodes[1], to = nodes[3]),
+        0.01,
+        0.12,
+        (from = 0.1, to = 0.1),
+        100,
+        1.04,
+    ),
+    PSY.Line(
+        "Line2",
+        true,
+        0.0,
+        0.0,
+        Arc(from = nodes[1], to = nodes[2]),
+        0.01,
+        0.12,
+        (from = 0.1, to = 0.1),
+        100,
+        1.04,
+    ),
+    PSY.Line(
+        "Line3",
+        true,
+        0.0,
+        0.0,
+        Arc(from = nodes[2], to = nodes[3]),
+        0.02,
+        0.9,
+        (from = 0.5, to = 0.5),
+        100,
+        1.04,
+    ),
+]
+
+branch_3bus_case9_fault(nodes) = [
+    PSY.Line(
+        "Line1",
+        true,
+        0.0,
+        0.0,
+        Arc(from = nodes[1], to = nodes[3]),
+        0.03,
+        0.36,
+        (from = 0.03, to = 0.03),
+        100,
+        1.04,
+    ),
+    PSY.Line(
+        "Line2",
+        true,
+        0.0,
+        0.0,
+        Arc(from = nodes[1], to = nodes[2]),
+        0.01,
+        0.12,
+        (from = 0.1, to = 0.1),
+        100,
+        1.04,
+    ),
+    PSY.Line(
+        "Line3",
+        true,
+        0.0,
+        0.0,
+        Arc(from = nodes[2], to = nodes[3]),
+        0.02,
+        0.9,
+        (from = 0.5, to = 0.5),
+        100,
+        1.04,
+    ),
+]
+
 ############### Load Data ########################
 
 loads_OMIB(nodes_OMIB) = [PSY.PowerLoad(
@@ -340,6 +418,12 @@ loads_3bus_case8(nodes_3bus) = [
     PSY.PowerLoad("Bus3", true, nodes_3bus[3], PSY.ConstantPower, 0.3, 0.1, 0.5, 0.3),
 ]
 
+loads_3bus_case9(nodes) = [
+    PowerLoad("Bus1", true, nodes[1], PowerSystems.ConstantPower, 0.5, 0.1, 1.5, 0.8),
+    PowerLoad("Bus2", true, nodes[2], PowerSystems.ConstantPower, 1.0, 0.3, 1.5, 0.8),
+    PowerLoad("Bus3", true, nodes[3], PowerSystems.ConstantPower, 0.3, 0.1, 0.5, 0.3),
+]
+
 ############### Infinite Sources Data ########################
 
 inf_gen_1_pu(nodes) = PSY.Source(
@@ -383,3 +467,20 @@ function get_admittance_matrix(nodes, branches)
     end
     return PSY.Ybus(sys2)[:, :]
 end
+
+function get_admittance_matrix_case9(nodes, branches)
+    sys2 = PSY.System(100.0, frequency = 60.0)
+    for bus in nodes
+        PSY.add_component!(sys2, bus)
+    end
+
+    #Add lines
+    for lines in branches
+        PSY.add_component!(sys2, lines)
+    end
+
+    #Make line 3 the dynamic line
+    make_dynamic_branch!(branches[3], sys2)
+
+    return PSY.Ybus(sys2)[:, :]
+end
diff --git a/test/old_tests/test_case9_dynbranches.jl b/test/old_tests/test_case9_dynbranches.jl
deleted file mode 100644
index eb2b105..0000000
--- a/test/old_tests/test_case9_dynbranches.jl
+++ /dev/null
@@ -1,310 +0,0 @@
-"""
-Case 9:
-This case study a three bus system with 1 machine (One d- One q-: 4th order model), a VSM of 19 states and an infinite source. The connection between buses 2 and 3 is modeled using dynamic lines.
-The perturbation trips two of the three circuits of line between buses 1 and 2, triplicating its impedance.
-"""
-
-##################################################
-############### LOAD DATA ########################
-##################################################
-
-############### Data Network ########################
-
-nodes_case9 = [
-    Bus(1, "Bus 1", "REF", 0, 1.02, (min = 0.94, max = 1.06), 138),
-    Bus(2, "Bus 2", "PV", 0, 1.00, (min = 0.94, max = 1.06), 138),
-    Bus(3, "Bus 3", "PQ", 0, 1.00, (min = 0.94, max = 1.06), 138),
-]
-
-branch_case9 = [
-    Line(
-        "Line1",
-        true,
-        0.0,
-        0.0,
-        Arc(from = nodes_case9[1], to = nodes_case9[3]),
-        0.01,
-        0.12,
-        (from = 0.1, to = 0.1),
-        100,
-        1.04,
-    ),
-    Line(
-        "Line2",
-        true,
-        0.0,
-        0.0,
-        Arc(from = nodes_case9[1], to = nodes_case9[2]),
-        0.01,
-        0.12,
-        (from = 0.1, to = 0.1),
-        100,
-        1.04,
-    ),
-]
-
-#Trip of Line 1.
-branch_case9_fault = [
-    Line(
-        "Line1",
-        true,
-        0.0,
-        0.0,
-        Arc(from = nodes_case9[1], to = nodes_case9[3]),
-        0.03,
-        0.36,
-        (from = 0.03, to = 0.03),
-        100,
-        1.04,
-    ),
-    Line(
-        "Line2",
-        true,
-        0.0,
-        0.0,
-        Arc(from = nodes_case9[1], to = nodes_case9[2]),
-        0.01,
-        0.12,
-        (from = 0.1, to = 0.1),
-        100,
-        1.04,
-    ),
-]
-
-loads_case9 = [
-    PowerLoad("Bus1", true, nodes_case9[1], PowerSystems.ConstantPower, 0.5, 0.1, 1.5, 0.8),
-    PowerLoad("Bus2", true, nodes_case9[2], PowerSystems.ConstantPower, 1.0, 0.3, 1.5, 0.8),
-    PowerLoad("Bus3", true, nodes_case9[3], PowerSystems.ConstantPower, 0.3, 0.1, 0.5, 0.3),
-]
-
-dyn_branch9 = [LITS.DynLine(
-    "Line3",
-    true,
-    Arc(from = nodes_case9[2], to = nodes_case9[3]),
-    0.02,
-    0.9,
-    (from = 0.5, to = 0.5),
-)]
-
-############### Data devices ########################
-
-inf_gen_case9 = StaticSource(
-    1, #number
-    :InfBus, #name
-    nodes_case9[1],#bus
-    1.00, #VR
-    0.0, #VI
-    0.000005,
-) #Xth
-
-######## Machine Data #########
-
-### Case 2: 4th Order Model with AVR (3-bus case) ###
-case9_machine = OneDOneQMachine(
-    0.0, #R
-    1.3125, #Xd
-    1.2578, #Xq
-    0.1813, #Xd_p
-    0.25, #Xq_p
-    5.89, #Td0_p
-    0.6, #Tq0_p
-    100.0,
-)   #MVABase
-
-######## Shaft Data #########
-
-### Shafts for Gen ###
-case9_shaft = SingleMass(
-    3.01, #H (M = 6.02 -> H = M/2)
-    0.0,
-) #D
-
-######## PSS Data #########
-cases_no_pss = PSSFixed(0.0)
-
-######## TG Data #########
-
-### No TG for Cases 1, 2, 3, 4 ###
-case9_no_tg = TGFixed(1.0) #eff
-
-########  AVR Data #########
-### AVRs for Case 2, 3, 4 and 5 ###
-case9_avr = AVRTypeI(
-    20.0, #Ka - Gain
-    0.01, #Ke
-    0.063, #Kf
-    0.2, #Ta
-    0.314, #Te
-    0.35, #Tf
-    0.001, #Tr
-    5.0, #Vrmax
-    -5.0, #Vrmin
-    0.0039, #Ae - 1st ceiling coefficient
-    1.555,
-) #Be - 2nd ceiling coefficient
-
-### Case 7 Generators ###
-case9_gen = PSY.DynamicGenerator(
-    1, #Number
-    :Case9Gen,
-    nodes_case9[2], #bus
-    1.0, # ω_ref,
-    1.0124, #V_ref
-    0.6, #P_ref
-    case9_machine, #machine
-    case9_shaft, #shaft
-    case9_avr, #avr
-    case9_no_tg, #tg
-    cases_no_pss,
-) #pss
-
-############### Inverter Data ########################
-
-converter = AvgCnvFixedDC(
-    138.0, #Rated Voltage
-    100.0,
-) #Rated MVA
-
-dc_source = FixedDCSource(1500.0) #Not in the original data, guessed.
-
-filt = LCLFilter(
-    0.08, #Series inductance lf in pu
-    0.003, #Series resitance rf in pu
-    0.074, #Shunt capacitance cf in pu
-    0.2, #Series reactance rg to grid connection (#Step up transformer or similar)
-    0.01,
-) #Series resistance lg to grid connection (#Step up transformer or similar)
-
-pll = PLL(
-    500.0, #ω_lp: Cut-off frequency for LowPass filter of PLL filter.
-    0.084, #k_p: PLL proportional gain
-    4.69,
-) #k_i: PLL integral gain
-
-virtual_H = VirtualInertia(
-    2.0, #Ta:: VSM inertia constant
-    400.0, #kd:: VSM damping coefficient
-    20.0, #kω:: Frequency droop gain in pu
-    2 * pi * 50.0,
-) #ωb:: Rated angular frequency
-
-Q_control = ReactivePowerDroop(
-    0.2, #kq:: Reactive power droop gain in pu
-    1000.0,
-) #ωf:: Reactive power cut-off low pass filter frequency
-
-outer_control = VirtualInertiaQdroop(virtual_H, Q_control)
-
-vsc = CombinedVIwithVZ(
-    0.59, #kpv:: Voltage controller proportional gain
-    736.0, #kiv:: Voltage controller integral gain
-    0.0, #kffv:: Binary variable enabling the voltage feed-forward in output of current controllers
-    0.0, #rv:: Virtual resistance in pu
-    0.2, #lv: Virtual inductance in pu
-    1.27, #kpc:: Current controller proportional gain
-    14.3, #kiv:: Current controller integral gain
-    0.0, #kffi:: Binary variable enabling the current feed-forward in output of current controllers
-    50.0, #ωad:: Active damping low pass filter cut-off frequency
-    0.2,
-) #kad:: Active damping gain
-
-case9_inv = PSY.DynamicInverter(
-    2, #number
-    :DARCO, #name
-    nodes_case9[3], #bus location
-    1.0, #ω_ref
-    0.8, #V_ref
-    0.5, #P_ref
-    -0.3, #Q_ref
-    100.0, #MVABase
-    converter, #Converter
-    outer_control, #OuterControl
-    vsc, #Voltage Source Controller
-    dc_source, #DC Source
-    pll, #Frequency Estimator
-    filt,
-) #Output Filter
-
-######################### Dynamical System ########################
-
-case9_DynSystem = PSY.System(
-    nodes_case9,
-    branch_case9,
-    [case9_inv, case9_gen],
-    vcat(inf_gen_case9, loads_case9),
-    100.0,
-    50.0,
-    dyn_branch9,
-);
-
-##################################################
-############### SOLVE PROBLEM ####################
-##################################################
-
-dx0 = zeros(LITS.get_total_rows(case9_DynSystem))
-x0 = [
-    1.00, #V1_R
-    1.00, #V2_R
-    1.00, #V3_R
-    0.0, #V1_I
-    -0.01, #V2_I
-    -0.01, #V3_I
-    0.0, #δω_vsm
-    0.2, #δθ_vsm
-    0.025, #qm
-    0.0015, #ξ_d
-    -0.07, #ξ_q
-    0.05, #γ_d
-    -0.001, #γ_q
-    0.95, #ϕ_d
-    -0.10, #ϕ_q
-    1.004, #vpll_d
-    0.0, #vpll_q
-    0.0, #ε_pll
-    0.1, #δθ_pll
-    0.5, #id_cv
-    0.0, #iq_cv
-    0.95, #vod
-    -0.1, #voq
-    0.49, #iod
-    -0.1, #ioq
-    1.0, #eq_p
-    0.47, #ed_p
-    0.6, #δ
-    1.0, #ω
-    2.1, #Vf
-    0.28, #Vr1
-    -0.39, #Vr2,
-    1.0, #Vm
-    0.5, #IL1_R
-    0.5,
-] #IL1_I
-
-case9_inv.inner_vars[13] = 0.95 #Vd_cnv var
-case9_inv.inner_vars[14] = -0.1 #Vq_cnv var
-Ybus_fault = Ybus(branch_case9_fault, nodes_case9)[:, :]
-u0 = [0.2]
-tspan = (0.0, 30.0);
-
-#Find initial condition
-inif! = (out, x) -> LITS.system!(out, dx0, x, (Ybus_fault, case9_DynSystem), 0.0)
-sys_solve = nlsolve(inif!, x0, xtol = :1e-8, ftol = :1e-8, method = :trust_region)
-x0_init = sys_solve.zero
-
-#Define Fault using Callbacks
-cb = DiffEqBase.DiscreteCallback(LITS.change_t_one, LITS.Y_change!)
-
-#Define Simulation Problem
-sim = Simulation(case9_DynSystem, tspan, Ybus_fault, cb, x0_init)
-
-#Solve problem in equilibrium
-run_simulation!(sim, IDA());
-
-#Obtain data for voltages
-series = get_voltagemag_series(sim, 2)
-zoom = [
-    (series[1][ix], series[2][ix])
-    for (ix, s) in enumerate(series[1]) if (s > 0.90 && s < 1.6)
-]
-
-@test sim.solution.retcode == :Success
diff --git a/test/test_case9_dynbranches.jl b/test/test_case9_dynbranches.jl
new file mode 100644
index 0000000..d52bfe6
--- /dev/null
+++ b/test/test_case9_dynbranches.jl
@@ -0,0 +1,121 @@
+"""
+Case 9:
+This case study a three bus system with 1 machine (One d- One q-: 4th order model), a VSM of 19 states and an infinite source. The connection between buses 2 and 3 is modeled using dynamic lines.
+The perturbation trips two of the three circuits of line between buses 1 and 2, triplicating its impedance.
+"""
+
+##################################################
+############### LOAD DATA ########################
+##################################################
+
+############### Data Network ########################
+
+nodes_case9 = nodes_3bus()
+
+branch_case9 = branch_3bus_case9(nodes_case9)
+
+#Trip of Line 1.
+branch_case9_fault = branch_3bus_case9_fault(nodes_case9)
+
+loads_case9 = loads_3bus_case9(nodes_case9)
+
+############### Data devices ########################
+
+inf_gen_case9 = inf_gen_1_pu(nodes_case9)
+
+######## Machine Data #########
+
+### Case 9 Generators ###
+case9_gen = dyn_gen_case9(nodes_case9)
+
+############### Inverter Data ########################
+
+case9_inv = inv_case78(nodes_case9)
+
+######################### Dynamical System ########################
+
+sys = system_case9(
+    nodes_case9,
+    branch_case9,
+    loads_case9,
+    [inf_gen_case9],
+    [case9_inv],
+    [case9_gen],
+);
+
+##################################################
+############### SOLVE PROBLEM ####################
+##################################################
+
+#Compute Y_bus after fault
+Ybus_fault = get_admittance_matrix_case9(nodes_case9, branch_case9_fault)
+
+#Time span
+tspan = (0.0, 30.0)
+
+#Initial guess
+x0_guess = [
+    1.00, #V1_R
+    1.00, #V2_R
+    1.00, #V3_R
+    0.0, #V1_I
+    -0.01, #V2_I
+    -0.01, #V3_I
+    0.0, #δω_vsm
+    0.2, #δθ_vsm
+    0.025, #qm
+    0.0015, #ξ_d
+    -0.07, #ξ_q
+    0.05, #γ_d
+    -0.001, #γ_q
+    0.95, #ϕ_d
+    -0.10, #ϕ_q
+    1.004, #vpll_d
+    0.0, #vpll_q
+    0.0, #ε_pll
+    0.1, #δθ_pll
+    0.5, #id_cv
+    0.0, #iq_cv
+    0.95, #vod
+    -0.1, #voq
+    0.49, #iod
+    -0.1, #ioq
+    1.0, #eq_p
+    0.47, #ed_p
+    0.6, #δ
+    1.0, #ω
+    2.1, #Vf
+    0.28, #Vr1
+    -0.39, #Vr2,
+    1.0, #Vm
+    0.5, #IL1_R
+    0.5,
+] #IL1_I
+
+#Define Fault: Change of YBus
+Ybus_change = LITS.ThreePhaseFault(
+    1.0, #change at t = 1.0
+    Ybus_fault,
+) #New YBus
+
+sim = Simulation(
+    sys, #system
+    tspan, #time span
+    Ybus_change, #Type of perturbation
+    initial_guess = x0_guess, #initial guess.
+)
+
+#Run simulation
+run_simulation!(
+    sim, #simulation structure
+    IDA(),
+) #Sundials DAE Solver
+
+#Obtain data for voltages
+series = get_voltagemag_series(sim, 2)
+zoom = [
+    (series[1][ix], series[2][ix])
+    for (ix, s) in enumerate(series[1]) if (s > 0.90 && s < 1.6)
+]
+
+@test sim.solution.retcode == :Success