diff --git a/README.md b/README.md
index 76574650..e4ad0ac7 100644
--- a/README.md
+++ b/README.md
@@ -11,7 +11,7 @@
 
 ## Version Advisory
 
-- PowerNetworkMatrices.jl will work with Julia v1.6+.
+- PowerNetworkMatrices.jl will work with Julia v1.6+. In Mac the requirement is to have Julia 1.9.2+ due to issues with MKL and dynamic library loading.
 - PowerNetworkMatrices.jl exports Matrix methods that were available in PowerSystems.jl version 1.0 and have been implemented as a separate package.
 
 ## License
diff --git a/docs/make.jl b/docs/make.jl
index f0497209..bce9bbbc 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -4,7 +4,13 @@ import DataStructures: OrderedDict
 pages = OrderedDict(
     "Welcome Page" => "index.md",
     "Quick Start Guide" => "quick_start_guide.md",
-    "Tutorials" => "tutorials/intro_page.md",
+    "Tutorials" => Any[
+        "Incidence, BA and ABA matrices" => "tutorials/tutorial_Incidence_BA_ABA_matrices.md",
+        "PTDF matrix" => "tutorials/tutorial_PTDF_matrix.md",
+        "VirtualPTDF matrix" => "tutorials/tutorial_VirtualPTDF_matrix.md",
+        "LODF matrix" => "tutorials/tutorial_LODF_matrix.md",
+        "VirtualLODF matrix" => "tutorials/tutorial_VirtualLODF_matrix.md",
+    ],
     "Public API Reference" => "api/public.md",
     "Internal API Reference" => "api/internal.md",
 )
@@ -15,7 +21,7 @@ makedocs(;
         mathengine = Documenter.MathJax(),
         prettyurls = haskey(ENV, "GITHUB_ACTIONS")),
     sitename = "PowerNetworkMatrices.jl",
-    authors = "Jose Daniel Lara, Alessandro Castelli, Sourabh Dalvi",
+    authors = "Jose Daniel Lara, Alessandro Francesco Castelli, Sourabh Dalvi",
     pages = Any[p for p in pages],
     clean = true,
 )
diff --git a/docs/src/assets/VirtualLODF_scheme.png b/docs/src/assets/VirtualLODF_scheme.png
new file mode 100644
index 00000000..7c15d66a
Binary files /dev/null and b/docs/src/assets/VirtualLODF_scheme.png differ
diff --git a/docs/src/assets/VirtualPTDF_scheme.png b/docs/src/assets/VirtualPTDF_scheme.png
new file mode 100644
index 00000000..41ac7640
Binary files /dev/null and b/docs/src/assets/VirtualPTDF_scheme.png differ
diff --git a/docs/src/index.md b/docs/src/index.md
index 220a673c..c4b3a227 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -6,7 +6,11 @@ CurrentModule = PowerNetworkMatrices
 
 ## Overview
 
-`PowerNetworkMatrices.jl` documentation and code are organized according to the needs of different
+`PowerNetworkMatrices.jl` is a [`Julia`](http://www.julialang.org) package for
+the evaluation of network matrices given the system's data. The package allows to compute 
+the matrices according to different methods, providing a flexibe and powerful tool.
+
+The documentation and code are organized according to the needs of different
 users depending on their skillset and requirements. In broad terms there are three categories:
 
 - **Modeler**: Users that want to run a particular analysis or experiment and use `PowerNetworkMatrices.jl` to develop data sets.
diff --git a/docs/src/quick_start_guide.md b/docs/src/quick_start_guide.md
index f2fc45f3..52ab2467 100644
--- a/docs/src/quick_start_guide.md
+++ b/docs/src/quick_start_guide.md
@@ -1 +1,42 @@
 # Quick Start Guide
+
+!!! note
+    `PowerSystemCaseBuilder.jl` is a helper library that makes it easier to reproduce examples in the documentation and tutorials. Normally you would pass your local files to create the system data instead of calling the function `build_system`.
+    For more details visit [PowerSystemCaseBuilder Documentation](https://nrel-sienna.github.io/PowerSystems.jl/stable/tutorials/powersystembuilder/)
+
+For more details about loading data and adding more dynamic components check the
+[Creating a System with Dynamic devices](https://nrel-sienna.github.io/PowerSystems.jl/stable/modeler_guide/system_dynamic_data/)
+section of the documentation in `PowerSystems.jl`.
+
+## Loading data
+
+Data can be loaded from a pss/e raw file and a pss/e dyr file.
+
+``` @repl quick_start_guide
+using PowerNetworkMatrices
+using PowerSystemCaseBuilder
+
+const PNM = PowerNetworkMatrices
+const PSB = PowerSystemCaseBuilder
+
+sys = PSB.build_system(PSB.PSITestSystems, "c_sys5")
+```
+
+## Computation of the PTDF matrix
+
+Once system data is loaded, netwrok matrices can be evaluated. The following
+example shows how the PTDF matrix is computed.
+
+The function `PTDF` is called for the evaluation of the matrix and other data. These
+are stored in a structure of type `PTDF`.
+
+``` @repl quick_start_guide
+# evaluate the PTDF structure containing the matrix and other data.
+ptdf_matrix = PNM.PTDF(sys);
+
+# show the PTDF matrix.
+PNM.get_data(ptdf_matrix)
+```
+
+As it can be seen, PTDF matrix is stored such that the number of rows is equal
+to the number of buses, number of columns equal to the number of branches.
diff --git a/docs/src/tutorials/tutorial_Incidence_BA_ABA_matrices.md b/docs/src/tutorials/tutorial_Incidence_BA_ABA_matrices.md
new file mode 100644
index 00000000..ee7c215e
--- /dev/null
+++ b/docs/src/tutorials/tutorial_Incidence_BA_ABA_matrices.md
@@ -0,0 +1,145 @@
+# Incidence, BA and ABA matrices
+
+In this tutorial the `IncidenceMatrix`, `BA_matrix` and `ABA_matrix` are presented. 
+The methods used for their evaluation, as well as how data is stored is shown in 
+the following subsections.
+
+The matrices here presented are the building blocks for the compuation of the PTDF and LODF matrices.
+
+## IncidenceMatrix
+
+The `PowerNetworkMatrices` package defines the structure `IncidenceMatrix`, which 
+store the Incidence Matrix of the considered system as well as the most relevant network data.
+
+At first, the `System` data is loaded.
+
+``` @repl tutorial_Incidence_BA_ABA_matrices
+using PowerNetworkMatrices
+using PowerSystemCaseBuilder
+
+const PNM = PowerNetworkMatrices
+const PSB = PowerSystemCaseBuilder
+
+sys = PSB.build_system(PSB.PSITestSystems, "c_sys5");
+```
+
+Then the Incidence Matrix is computed as follows:
+
+``` @repl tutorial_Incidence_BA_ABA_matrices
+incidence_matrix = PNM.IncidenceMatrix(sys);
+```
+
+The `incidence_matrix` variable is a structure of type `IncidenceMatrix` 
+featuring the following fields:
+
+``` @repl tutorial_Incidence_BA_ABA_matrices
+# axis names: row and column names.
+# row names: names of the branches
+# column names: names of the buses
+incidence_matrix.axes
+
+# data: Incidence Matrix
+incidence_matrix.data
+
+# lookup: dictionary linking the branche names and bus numbers with the row 
+# and column numbers, respectively.
+incidence_matrix.axes
+
+# ref_bus_positions: set containing the positions of the reference buses.
+# this represents the positions where to add the column of zeros. Please refer to the 
+# exaple in the BA matrix for more details.
+incidence_matrix.ref_bus_positions
+```
+
+Please note that the matrix data can be easily access by using the following function:
+
+``` @repl tutorial_Incidence_BA_ABA_matrices
+PNM.get_data(incidence_matrix)
+```
+
+Note that the number of columns is lower than the actual number of system buses since 
+the column related to the reference bus is discarded.
+
+## BA_Matrix
+
+The `BA_Matrix` is a structure containing the matrix coming from the product of the
+`IncidenceMatrix` and the diagonal matrix contianing the impedence of the system's branches ("B" matrix).
+
+The `BA_Matrix` is computed as follows:
+
+``` @repl tutorial_Incidence_BA_ABA_matrices
+ba_matrix = PNM.BA_Matrix(sys);
+```
+
+As for the `IncidenceMatrix`, here too the `BA_Matrix` structure feature the same fields.
+
+An example related to accessing the matrix data is now provided:
+``` @repl tutorial_Incidence_BA_ABA_matrices
+# access data by explicitly calling the field "data"
+ba_matrix.data
+
+# or by using the "get_data" function
+PNM.get_data(ba_matrix)
+```
+Note that the number of columns is lower than the actual number of system buses since 
+the column related to the reference bus is discarded.
+
+To add the column of zeros related to the reference bus, it is necessary to use the 
+information contained in the `ref_bus_positions` field.
+
+``` @repl tutorial_Incidence_BA_ABA_matrices
+new_ba_matrix = hcat(
+    ba_matrix.data[:,1:collect(ba_matrix.ref_bus_positions)[1]-1], 
+    zeros(size(ba_matrix, 1), 1), 
+    ba_matrix.data[:, collect(ba_matrix.ref_bus_positions)[1]:end]
+    )
+```
+
+However, trying to change the data field with a matrix of different dimension 
+will result in an error.
+
+``` @repl tutorial_Incidence_BA_ABA_matrices
+ba_matrix.data = hcat(
+    ba_matrix.data[:,1:collect(ba_matrix.ref_bus_positions)[1]-1], 
+    zeros(size(ba_matrix, 1), 1), 
+    ba_matrix.data[:, collect(ba_matrix.ref_bus_positions)[1]:end]
+    )
+```
+
+
+## ABA_Matrix
+
+The `ABA_Matrix` is a structure containing the matrix coming from the product of the
+`IncidenceMatrix` and the `BA_Matrix`.
+It features the same fields as the `IncidenceMatrix` and the `BA_Matrix`, plus the `K` one.
+The field `ABA_Matrix.K` stores the LU factorization matrices (using the 
+methods contained in the package `KLU`).
+
+To evaluate the `ABA_Matrix`, the following command is sufficient:
+
+``` @repl tutorial_Incidence_BA_ABA_matrices
+aba_matrix = ABA_Matrix(sys);
+```
+
+By default the LU factorization matrices are not computed, leaving the `K` field empty.
+In case these are wanted, the keyword `factorize` must be true.
+
+``` @repl tutorial_Incidence_BA_ABA_matrices
+aba_matrix_with_LU = ABA_Matrix(sys; factorize=true);
+
+aba_matrix_with_LU.K
+```
+
+If the `ABA_Matrix` is already computed but the LU factorization was not performed, this can be done by considering the following command:
+
+``` @repl tutorial_Incidence_BA_ABA_matrices
+aba_matrix.K
+aba_matrix = factorize(aba_matrix);
+aba_matrix.K
+```
+
+The following command can then be used to check if the `ABA_Matrix` contains the LU factorization matrices:
+
+``` @repl tutorial_Incidence_BA_ABA_matrices
+is_factorized(aba_matrix_new)
+```
\ No newline at end of file
diff --git a/docs/src/tutorials/tutorial_LODF_matrix.md b/docs/src/tutorials/tutorial_LODF_matrix.md
new file mode 100644
index 00000000..348edc35
--- /dev/null
+++ b/docs/src/tutorials/tutorial_LODF_matrix.md
@@ -0,0 +1,100 @@
+# LODF matrix
+
+In this tutorial the methods for computing the Line Outage Distribution Factor (`LODF`) are presented.
+Before diving into this tutorial we encourage the user to load `PowerNetworkMatrices`, hit the `?` key in the REPL terminal and look for the documentiont of the different `LODF` methods avialable.
+
+## Evaluation of the `LODF` matrix
+
+As for the `PTDF` matrix, the `LODF` one can be evaluated according to two different approaches:
+- `Dense`: considers functions for dense matrix multiplication and inversion
+- `KLU`: considers functions for sparse matrix multiplication  and inversion(**default**)
+
+The evaluation of the `LODF` matrix can be easily performed starting from importing the system's data and then by simply calling the `LODF` method.
+
+``` @repl tutorial_PTDF_matrix
+using PowerNetworkMatrices
+using PowerSystemCaseBuilder
+
+const PNM = PowerNetworkMatrices
+const PSB = PowerSystemCaseBuilder
+
+# get the System data
+sys = PSB.build_system(PSB.PSITestSystems, "c_sys5");
+
+# compute the LODF matrix
+lodf_1 = LODF(sys);
+
+lodf_2 = LODF(sys, linear_solver="Dense");
+
+# show matrix
+get_lodf_data(lodf_1)
+```
+
+Advanced users might be interested in computing the `LODF` matrix starting from either the `branches` and `buses` data (`CASE 1`), the `IncidenceMatrix` and `PTDF` structures (`CASE 2`), or by the information related to `IncidenceMatrix`, `BA_Matrix` and `ABA_Matrix` (`CASE 3`).
+
+``` @repl tutorial_PTDF_matrix
+# CASE 1
+
+# get the branches and buses data
+branches = PNM.get_ac_branches(sys);
+buses = PNM.get_buses(sys);
+
+# compute the LODF matrix from branches and buses data
+lodf_3 = LODF(branches, buses);
+
+# CASE 2
+
+# get the Incidence and PTDF matrix
+a = IncidenceMatrix(sys);
+ptdf = PTDF(sys);
+
+# compute LODF matrix with the two netwrok matrices
+lodf_4 = LODF(a, ptdf);
+
+# CASE 3
+
+# get the BA and ABA matrices (ABA matrix must include LU factorization 
+# matrices)
+ba = BA_Matrix(sys);
+aba = ABA_Matrix(sys, factorize = true);
+
+# compute LODF matrix with the three netwrok matrices
+lodf_5 = LODF(a, aba, ba);
+```
+
+**NOTE:** whenever the method `LODF(sys::System)` is used, the methods previously defined for `CASE 1` and `CASE 2` are executed in sequence. Therefore the method `LODF(a::IncidenceMatrix, ptdf::PTDF)` is the default one when evaluating the `LODF` matrix from the `System` data directly.
+
+
+## Available methods for the computation of the `LODF` matrix
+
+For those methods that either require the evaluation of the `PTDF` matrix, or that execute this evaluation internally, two different approaches casen be used.
+
+As for the `PTDF` matrix, here too the optional argument `linear_solver` can be specified with either `KLU` (for spars matrix calculation) or `Dense` (for sparse matrix calculation).
+
+``` @repl tutorial_PTDF_matrix
+lodf_dense = LODF(sys, linear_solver="Dense");
+```
+
+**NOTE (1):** by default the "KLU" method is selected, which appeared to require significant less time and memory with respect to "Dense".
+Please note that wether the `KLU` or `Dense` method is used, the resultig `LODF` matrix is stored as a dense one.
+
+**NOTE (2):** for the moment, the method `LODF(a::IncidenceMatrix, aba::ABA_Matrix, ba::BA_Matrix)` will take `KLU` as `linear_solver` option.
+
+## "Sparse" `LODF` matrix
+
+The `LODF` matrix can be computed in a "sparse" fashion by defining the input argument `tol`. If this argument is defined, then elements of the `LODF` matrix whose absolute values are below the set threshold are dropped. In addition, the matrix will be stored as a sparse one of type `SparseArrays.SparseMatrixCSC{Float64, Int64}` type instead of `Matrix{Float64}` one.
+
+By considering an "extreme" value of 0.4 as `tol`, the `LODF` matrix can be computed as follows:
+
+``` @repl tutorial_PTDF_matrix
+lodf_sparse = LODF(sys, tol=0.4);
+get_lodf_data(lodf_sparse)
+```
+
+Please consider that 0.4 was used for the purpose of this tutorial. In practice much smaller values are used (e.g., 1e-5).
+
+**NOTE (1):** elements whose absolute values exceed the `tol` argument are removed from the `LODF` matrix *after* this has been computed.
+
+**NOTE (2):** the `tol` argument does not refer to the "sparsification" tolerance of the `PTDF` matrix that is computed in the `LODF` method.
+
+**NOTE (3):** in case the method `LODF(a::IncidenceMatrix, ptdf::PTDF)` is considerd, an error will be thrown whenever the `tol` argument in the `PTDF` structure used as input is different then 1e-15.
\ No newline at end of file
diff --git a/docs/src/tutorials/tutorial_PTDF_matrix.md b/docs/src/tutorials/tutorial_PTDF_matrix.md
new file mode 100644
index 00000000..ac286691
--- /dev/null
+++ b/docs/src/tutorials/tutorial_PTDF_matrix.md
@@ -0,0 +1,101 @@
+# PTDF matrix
+
+In this tutorial the methods for computing the Power Transfer Distribution Factors (`PTDF`) are presented.
+Before diving into this tutorial we encourage the user to load `PowerNetworkMatrices`, hit the `?` key in the REPL terminal and look for the documentiont of the different `PTDF` methods avialable.
+
+## Evaluation of the `PTDF` matrix
+
+The `PTDF` matrix can be evaluated according to two different approaches:
+- `Dense`: considers functions for dense matrix multiplication and inversion
+- `KLU`: considers functions for sparse matrix multiplication and inversion (**default**)
+
+The evaluation of the `PTDF` matrix can be easily performed starting from importing the system's data and then by simply calling the `PTDF` method.
+
+``` @repl tutorial_PTDF_matrix
+using PowerNetworkMatrices
+using PowerSystemCaseBuilder
+
+const PNM = PowerNetworkMatrices
+const PSB = PowerSystemCaseBuilder
+
+sys = PSB.build_system(PSB.PSITestSystems, "c_sys5");
+
+ptdf_1 = PTDF(sys);
+
+get_ptdf_data(ptdf_1)
+```
+
+Advanced users might be interested in computing the `PTDF` matrix starting from either the data contained in the `IncidenceMatrix` and `BA_matrix` structures, or by the information related to the `branches` and `buses` of the system.
+
+``` @repl tutorial_PTDF_matrix
+# evaluate the BA_matrix and Incidence_Matrix
+ba_matrix = BA_Matrix(sys);
+a_matrix = IncidenceMatrix(sys);
+
+# get the PTDF matrix starting from the values of the 
+# previosly cumputed matrices
+ptdf_2 = PTDF(a_matrix, ba_matrix);
+get_ptdf_data(ptdf_2)
+
+# get the buses and branches of the system
+branches = PNM.get_ac_branches(sys);
+buses = PNM.get_buses(sys);
+ptdf_3 = PTDF(branches, buses);
+get_ptdf_data(ptdf_3)
+```
+
+NOTE: both the `get_ac_branches` and `get_ac_branches` functions are not exported by the `PowerNetworkMatrices` package, and therefore require the package name to be called as a prefix. However, they are shown here just for the sake of making an example.
+
+## Available methods for the computation of the `PTDF` matrix
+
+As previously mentioned, the `PTDF` matrix can be evaluated considering different approaches. The method can be selected by specifying the field `linear_solver` in the `PTDF` function.
+
+``` @repl tutorial_PTDF_matrix
+ptdf_dense = PTDF(sys, linear_solver="Dense");
+get_ptdf_data(ptdf_dense)
+
+ptdf_klu = PTDF(sys, linear_solver="KLU");
+get_ptdf_data(ptdf_klu)
+```
+
+By default the "KLU" method is selected, which appeared to require significant less time and memory with respect to "Dense".
+Please note that either the `KLU` or `Dense` method isi used, the resultig `PTDF` matrix is stored as a dense one.
+
+## Evaluating the `PTDF` matrix considering distributed slack bus
+
+Whenever needed, the `PTDF` matrix can be computed with a distributed slack bus. To do so, a vecotr of type `Vector{Float64}` needs to be defined, specifying the weights for each bus of the system. These weights identify how the load on the slakc bus is redistributed accross the system.
+
+``` @repl tutorial_PTDF_matrix
+# consider equal distribution accross each bus for this example
+buscount = length(PNM.get_buses(sys));
+dist_slack = 1 / buscount * ones(buscount);
+dist_slack_array = dist_slack / sum(dist_slack);
+```
+
+Once the vector of the weights is defined, the `PTDF` matrix can be computed by defining the input argument `dist_slack` (empty array `Float64[]` by default):
+
+``` @repl tutorial_PTDF_matrix
+ptdf_distr = PTDF(sys, dist_slack=dist_slack_array);
+```
+
+The diffrence between a the matrix computed with and without the `dist_slack` field defined can be seen as follows:
+
+``` @repl tutorial_PTDF_matrix
+# with no distributed slack bus
+get_ptdf_data(ptdf_klu)
+# with distributed slack bus
+get_ptdf_data(ptdf_distr)
+```
+
+## "Sparse" `PTDF` matrix
+
+The `PTFD` matrix can be computed in a "sparse" fashion by defining the input argument `tol`. If this argument is defined, then elements of the `PTDF` matrix whose absolute values are below the set threshold are dropped. In addition, the matrix will be stored as a sparse one of type `SparseArrays.SparseMatrixCSC{Float64, Int64}` instead of `Matrix{Float64}`.
+
+By considering an "extreme" value of 0.2 as `tol`, the `PTDF` matrix can be computed as follows:
+
+``` @repl tutorial_PTDF_matrix
+ptdf_sparse = PTDF(sys, tol=0.2);
+get_ptdf_data(ptdf_sparse)
+```
+
+**NOTE:** 0.2 was used for the purpose of this tutorial. In practice much smaller values are used (e.g., 1e-5).
\ No newline at end of file
diff --git a/docs/src/tutorials/tutorial_VirtualLODF_matrix.md b/docs/src/tutorials/tutorial_VirtualLODF_matrix.md
new file mode 100644
index 00000000..803feddc
--- /dev/null
+++ b/docs/src/tutorials/tutorial_VirtualLODF_matrix.md
@@ -0,0 +1,98 @@
+# VirtualLODF
+
+The `VirtualLODF` structure follows the same philosofy as the `VirtualPTDF`: it contains rows of the original `LODF` matrix, evaluated and cached on demand.
+
+Refer to the different arguments of the `VirtualLODF` methods by looking at the "Public API Reference" page.
+
+## How the `VirtualLODF` works
+
+The `VirtualLODF` structure retains many of the similarities of the `VirtualPTDF`. However, its computation is more complex and requires some additional data.
+
+Starting from the system data, the `IncidenceMatrix`, `BA_Matrix` and `ABA_Matrix` (with relative LU factorization matrices) are evaluated. The `ABA_Matrix` and `BA_Matrix` are used for the computation of the diagonal elements of the `PTDF` matrix, and this vector is stored in the `VirtualLODF` structure together with the other structures abovementioned.
+
+Once the `VirtualLODF` is initialized, each row of the matrix can be evaluated separately and on user request. The algorithmic procedure is the following:
+1. Define the `VirtualPTDF` structure
+2. Call any element of the matrix to define and store the relative row as well as showing the selected element
+
+Regarding point 2, if the row has been stored previosly then the desired element is just loaded from the cache and shown.
+
+The flowchart below shows how the `VirtualLODF` is structured and how it works. Examples will be presented in the following sections.
+
+```@raw html
+<img src="../assets/VirtualLODF_scheme.png"/>
+```
+
+## Initialize `VirtualLODF` and compute/access row/element
+
+As for the `LODF` matrix, at first the `System` data must be loaded. The "RTS GMLC" systems is considered as example:
+
+``` @repl tutorial_VirtualPTDF_matrix
+using PowerNetworkMatrices
+using PowerSystemCaseBuilder
+
+const PNM = PowerNetworkMatrices;
+const PSB = PowerSystemCaseBuilder;
+
+sys = PSB.build_system(PSB.PSISystems, "RTS_GMLC_DA_sys");
+```
+
+At this point the `VirtualLODF` is initialized with the following simple command:
+
+``` @repl tutorial_VirtualPTDF_matrix
+v_lodf = VirtualLODF(sys);
+```
+
+Now, an element of the matrix can be computed by calling the branch name and bus number:
+
+``` @repl tutorial_VirtualPTDF_matrix
+el_C31_2_105 = v_lodf["C31-2", "A3"]
+```
+
+This element represent the portion flowing on line "A3" now diverted on line "C31-2" as a consequence of its outage.
+
+Alternatively, the number of the branch and bus (corresponding to the number of the PTDF row and column) can be used. In this case the row and column numbers are mapped by the dictonaries contained in the `lookup` field. 
+
+``` @repl tutorial_VirtualPTDF_matrix
+row_number = v_lodf.lookup[1]["C31-2"]
+col_number = v_lodf.lookup[2]["A3"]
+el_C31_2_105_bis = v_lodf[row_number, col_number]
+```
+
+**NOTE**: this example was made for the sake of completeness and considering the actual branch names is reccomended.
+
+As previosly mentioned, in order to evalute a single element of the `VirtualLODF`, the entire row related to the selected branch must be considered. For this reason it is cached for later calls.
+This is evident by looking at the following example:
+
+``` @repl tutorial_VirtualPTDF_matrix
+sys_2k = PSB.build_system(PSB.PSYTestSystems, "tamu_ACTIVSg2000_sys");
+
+v_lodf_2k = VirtualLODF(sys_2k);
+
+# evaluate PTDF row related to branch "ODESSA 2 0  -1001-ODESSA 3 0  -1064-i_1"
+@time v_lodf_2k["ODESSA 2 0  -1001-ODESSA 3 0  -1064-i_1", 
+                "BRYAN 1 1   -8159-BRYAN 1 0   -8158-i_1"]
+
+# call same element after the row has been stored
+@time v_lodf_2k["ODESSA 2 0  -1001-ODESSA 3 0  -1064-i_1", 
+                "BRYAN 1 1   -8159-BRYAN 1 0   -8158-i_1"]
+```
+
+## "Sparse" `VirtualPTDF`
+
+Sparsification of each row can be achieved in the same fashion as for the `LODF` matrix, by removing those elements whose absolute values is below a certain threshold.
+
+As for the example show for the `LODF` matrix, here to a very high values of 0.4 is considered for the `tol` field. Again, this value is considered just for the sake of this example.
+
+``` @repl tutorial_VirtualPTDF_matrix
+# smaller system for the next examples
+sys_2 = PSB.build_system(PSB.PSITestSystems, "c_sys5");
+
+v_lodf_dense = VirtualLODF(sys_2);
+v_lodf_sparse = VirtualLODF(sys_2, tol=0.4);
+```
+
+Let's now evaluate the same row as before and compare the results:
+``` @repl tutorial_VirtualPTDF_matrix
+original_row = [v_lodf_dense["1", j] for j in v_lodf_dense.axes[2]]
+sparse_row = [v_lodf_sparse["1", j] for j in v_lodf_sparse.axes[2]]
+```
\ No newline at end of file
diff --git a/docs/src/tutorials/tutorial_VirtualPTDF_matrix.md b/docs/src/tutorials/tutorial_VirtualPTDF_matrix.md
new file mode 100644
index 00000000..7afcb5d6
--- /dev/null
+++ b/docs/src/tutorials/tutorial_VirtualPTDF_matrix.md
@@ -0,0 +1,119 @@
+# VirtualPTDF
+
+Contrary to the traditional `PTDF` matrix, the `VirtualPTDF` is a stucture contatining rows of the original matrix, related to specific system branches.
+The different rows of the `PTDF` matrix are cached in the `VirtualPTDF` structure as they are evaluated. This allows to keep just the portion of the original matrix which is of interest to the user, avoiding the unecessary computation of the whole matrix.
+
+Refer to the different arguments of the `VirtualPTDF` methods by looking at the "Public API Reference" page.
+
+## How the `VirtualPTDF` works
+
+The `VirtualPTDF` is a structure containing everything needed to compute any row of the PTDF matrix and store it. To do so, the `VirtualPTDF` must store the BA matrix (coming from the `BA_Matrix` struct) and the inverse of the ABA matrix (coming from `ABA_MAtrix` struct). In particular, `KLU` is used to get the LU factorization matrices of the ABA matrix and these ones are stored, avoid the inversion.
+
+Once the `VirtualPTDF` is initialized, each row of the PTDF matrix can be evaluated separately. The algorithmic procedure is the following:
+1. Define the `VirtualPTDF` structure
+2. Call any element of the matrix to define and store the relative row as well as showing the selected element
+
+Regarding point 2, if the row has been stored previosly then the desired element is just loaded from the cache and shown.
+
+The flowchart below shows how the `VirtualPTDF` is structured and how it works. Examples will be presented in the following sections.
+
+```@raw html
+<img src="../assets/VirtualPTDF_scheme.png"/>
+```
+
+## Initialize `VirtualPTDF` and compute/access row/element
+
+As for the `PTDF` matrix, at first the `System` data must be loaded. The "RTS GMLC" systems is considered as example:
+
+``` @repl tutorial_VirtualPTDF_matrix
+using PowerNetworkMatrices
+using PowerSystemCaseBuilder
+
+const PNM = PowerNetworkMatrices;
+const PSB = PowerSystemCaseBuilder;
+
+sys = PSB.build_system(PSB.PSISystems, "RTS_GMLC_DA_sys");
+```
+
+At this point the `VirtualPTDF` is initialized with the following simple command:
+
+``` @repl tutorial_VirtualPTDF_matrix
+v_ptdf = VirtualPTDF(sys);
+```
+
+Now, an element of the matrix can be computed by calling the branch name and bus number:
+
+``` @repl tutorial_VirtualPTDF_matrix
+el_C31_2_105 = v_ptdf["C31-2", 105]
+```
+
+Alternatively, the number of the branch and bus (corresponding to the number of the PTDF row and column) can be used. In this case the row and column numbers are mapped by the dictonaries contained in the `lookup` field. 
+
+``` @repl tutorial_VirtualPTDF_matrix
+row_number = v_ptdf.lookup[1]["C31-2"]
+col_number = v_ptdf.lookup[2][105]
+el_C31_2_105_bis = v_ptdf[row_number, col_number]
+```
+
+**NOTE**: this example was made for the sake of completeness and considering the actual branch name and bus number is reccomended.
+
+As previosly mentioned, in order to evalute a single element of the `VirtualPTDF`, the entire row related to the selected branch must be considered. For this reason it is cached in the `VirtualPTDF` structure for later calls.
+This is evident by looking at the following example:
+
+``` @repl tutorial_VirtualPTDF_matrix
+sys_2k = PSB.build_system(PSB.PSYTestSystems, "tamu_ACTIVSg2000_sys");
+
+v_ptdf_2k = VirtualPTDF(sys_2k);
+
+# evaluate PTDF row related to branch "ODESSA 2 0  -1001-ODESSA 3 0  -1064-i_1"
+@time v_ptdf_2k["ODESSA 2 0  -1001-ODESSA 3 0  -1064-i_1", 8155]
+
+# call same element after the row has been stored
+@time v_ptdf_2k["ODESSA 2 0  -1001-ODESSA 3 0  -1064-i_1", 8155]
+```
+
+## `VirtualPTDF` with distributed slack bus
+
+As for the `PTDF` matrix, here too each row can be evaluated considering distibuted slack buses.
+A vecotr of type `Vector{Float64}` is defined, specifying the weights for each bus of the system. 
+
+``` @repl tutorial_VirtualPTDF_matrix
+# smaller system for the next examples
+sys_2 = PSB.build_system(PSB.PSITestSystems, "c_sys5");
+
+# consider equal distribution accross each bus for this example
+buscount = length(PNM.get_buses(sys_2));
+dist_slack = 1 / buscount * ones(buscount);
+dist_slack_array = dist_slack / sum(dist_slack);
+```
+
+Now initialize the `VirtualPTDF` by defining the `dist_slack` field with the vector of weights previosly computed:
+
+``` @repl tutorial_VirtualPTDF_matrix
+v_ptdf_distr = VirtualPTDF(sys_2, dist_slack=dist_slack_array);
+v_ptdf_orig = VirtualPTDF(sys_2);
+```
+
+Now check the difference with the same row related to the branch `"1"` evaluated withou considering distributed slack bus.
+
+``` @repl tutorial_VirtualPTDF_matrix
+row_distr = [v_ptdf_distr["1", j] for j in v_ptdf_distr.axes[2]]
+row_original = [v_ptdf_orig["1", j] for j in v_ptdf_orig.axes[2]]
+```
+
+## "Sparse" `VirtualPTDF`
+
+Sparsification of each row can be achieved in the same fashion as for the `PTDF` matrix, by removing those elements whose absolute values is below a certain threshold.
+
+As for the example show for the `PTDF` matrix, here to a very high values of 0.2 is considered for the `tol` field. Again, this value is considered just for the sake of this example.
+
+``` @repl tutorial_VirtualPTDF_matrix
+v_ptdf_dense = VirtualPTDF(sys_2);
+v_ptdf_sparse = VirtualPTDF(sys_2, tol=0.2);
+```
+
+Let's now evaluate the same row as before and compare the results:
+``` @repl tutorial_VirtualPTDF_matrix
+original_row = [v_ptdf_dense["1", j] for j in v_ptdf_dense.axes[2]]
+sparse_row = [v_ptdf_sparse["1", j] for j in v_ptdf_sparse.axes[2]]
+```
\ No newline at end of file
diff --git a/src/virtual_lodf_calculations.jl b/src/virtual_lodf_calculations.jl
index e4f6e5fd..41cadc91 100644
--- a/src/virtual_lodf_calculations.jl
+++ b/src/virtual_lodf_calculations.jl
@@ -203,7 +203,7 @@ function _getindex(
 )
     # check if value is in the cache
     if haskey(vlodf.cache, row)
-        return vlodf.cache[row][column]
+        return vlodf.cache.temp_cache[row][column]
     else
 
         # evaluate the value for the LODF column
diff --git a/src/virtual_ptdf_calculations.jl b/src/virtual_ptdf_calculations.jl
index bf5e30f1..13deeee3 100644
--- a/src/virtual_ptdf_calculations.jl
+++ b/src/virtual_ptdf_calculations.jl
@@ -191,7 +191,7 @@ function _getindex(
 )
     # check if value is in the cache
     if haskey(vptdf.cache, row)
-        return vptdf.cache[row][column]
+        return vptdf.cache.temp_cache[row][column]
     else
         # evaluate the value for the PTDF column
         # Needs improvement
diff --git a/test/Project.toml b/test/Project.toml
index 4fc37933..0fe90594 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,9 +1,11 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
 PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd"
 PowerSystemCaseBuilder = "f00506e0-b84f-492a-93c2-c0a9afc4364e"
 PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"