diff --git a/.gitignore b/.gitignore
index 59c8955c..fa98f6cd 100644
--- a/.gitignore
+++ b/.gitignore
@@ -150,3 +150,29 @@ dmypy.json
 
 # Pyre type checker
 .pyre/
+
+# LaTex files
+## Do not ignore the documentation document
+!documentation/main.pdf
+
+## Core latex/pdflatex auxiliary files:
+*.aux
+*.lof
+*.log
+*.lot
+*.fls
+*.out
+*.toc
+*.fmt
+*.fot
+*.cb
+*.cb2
+.*.lb
+
+## Build tool auxiliary files:
+*.fdb_latexmk
+*.synctex
+*.synctex(busy)
+*.synctex.gz
+*.synctex.gz(busy)
+*.pdfsync
diff --git a/conmech/dynamics/dynamics.py b/conmech/dynamics/dynamics.py
index 8e1325c1..a479c8eb 100644
--- a/conmech/dynamics/dynamics.py
+++ b/conmech/dynamics/dynamics.py
@@ -89,7 +89,8 @@ def __init__(
         self.acceleration_operator: np.ndarray
         self.elasticity: np.ndarray
         self.viscosity: np.ndarray
-        self._w_matrix = None
+        self._w_matrix: np.ndarray = None
+        self._local_stifness_matrices: np.ndarray = None
         self.__relaxation: Optional[np.ndarray] = None
         self.__relaxation_tensor: Optional[float] = None
         self.thermal_expansion: np.ndarray
@@ -107,14 +108,19 @@ def remesh(self, boundaries_description, create_in_subprocess):
         super().mesh.remesh(boundaries_description, create_in_subprocess)
         self.reinitialize_matrices()
 
-    def reinitialize_matrices(self):
+    def reinitialize_matrices(self, elements_density: Optional[np.ndarray] = None):
         (
             self.element_initial_volume,
             self.volume_at_nodes,
             U,
             V,
             self._w_matrix,
+            self._local_stifness_matrices,
         ) = get_basic_matrices(elements=self.mesh.elements, nodes=self.moved_nodes)
+
+        if elements_density is not None:
+            self._w_matrix = self.asembly_w_matrix_with_density(elements_density)
+
         (
             self.acceleration_operator,
             self.elasticity,
@@ -165,6 +171,17 @@ def reinitialize_matrices(self):
                     free_indices=self.mesh.free_indices,
                 )
 
+    def asembly_w_matrix_with_density(self, elements_density: np.ndarray):
+        w_matrix = np.zeros_like(self._w_matrix)
+        for element_index, element in enumerate(self.mesh.elements):
+            for i, global_i in enumerate(element):
+                for j, global_j in enumerate(element):
+                    w_matrix[:, :, global_i, global_j] += (
+                        elements_density[element_index]
+                        * self._local_stifness_matrices[:, :, element_index, i, j]
+                    )
+        return w_matrix
+
     @property
     def with_temperature(self):
         return isinstance(self.body_prop, TemperatureBodyProperties)
diff --git a/conmech/dynamics/factory/_dynamics_factory_2d.py b/conmech/dynamics/factory/_dynamics_factory_2d.py
index 00218934..c967a19e 100644
--- a/conmech/dynamics/factory/_dynamics_factory_2d.py
+++ b/conmech/dynamics/factory/_dynamics_factory_2d.py
@@ -23,6 +23,11 @@ def get_edges_features_matrix_numba(elements, nodes):
         (FEATURE_MATRIX_COUNT, nodes_count, nodes_count), dtype=np.double
     )
     element_initial_volume = np.zeros(elements_count)
+    # Local stifness matrices (w[0, 0], w[0, 1], w[1, 0], w[1, 1]) per mesh element
+    # Detailed description can be found in [LSM] Local stifness matrix
+    local_stifness_matrices = np.empty(
+        (DIMENSION, DIMENSION, elements_count, element_size, element_size)
+    )
 
     for element_index in range(elements_count):  # TODO: #65 prange?
         element = elements[element_index]
@@ -40,9 +45,21 @@ def get_edges_features_matrix_numba(elements, nodes):
         )
         int_matrix = np.array(
             [
-                [0, 1 / 6 * element_nodes[:, 0].sum(), 1 / 6 * element_nodes[:, 1].sum()],
-                [1 / 6 * element_nodes[:, 0].sum(), 1 / 12 * int_sqr[0], 1 / 24 * int_xy],
-                [1 / 6 * element_nodes[:, 1].sum(), 1 / 24 * int_xy, 1 / 12 * int_sqr[1]],
+                [
+                    0,
+                    1 / 6 * element_nodes[:, 0].sum(),
+                    1 / 6 * element_nodes[:, 1].sum(),
+                ],
+                [
+                    1 / 6 * element_nodes[:, 0].sum(),
+                    1 / 12 * int_sqr[0],
+                    1 / 24 * int_xy,
+                ],
+                [
+                    1 / 6 * element_nodes[:, 1].sum(),
+                    1 / 24 * int_xy,
+                    1 / 12 * int_sqr[1],
+                ],
             ]
         )
         jacobian = np.abs(np.linalg.det(element_nodes[1:] - element_nodes[0]))
@@ -85,6 +102,8 @@ def get_edges_features_matrix_numba(elements, nodes):
 
                 w = [[i_d_phi * j_d_phi for j_d_phi in j_d_phi_vec] for i_d_phi in i_d_phi_vec]
 
+                local_stifness_matrices[:, :, element_index, i, j] = element_volume * np.asarray(w)
+
                 edges_features_matrix[:, element[i], element[j]] += element_volume * np.array(
                     [
                         volume_at_nodes,
@@ -100,7 +119,7 @@ def get_edges_features_matrix_numba(elements, nodes):
                 )
 
     # Performance TIP: we need only sparse, triangular matrix (?)
-    return edges_features_matrix, element_initial_volume
+    return edges_features_matrix, element_initial_volume, local_stifness_matrices
 
 
 @numba.njit
diff --git a/conmech/dynamics/factory/_dynamics_factory_3d.py b/conmech/dynamics/factory/_dynamics_factory_3d.py
index 1141230a..c6ef16c9 100644
--- a/conmech/dynamics/factory/_dynamics_factory_3d.py
+++ b/conmech/dynamics/factory/_dynamics_factory_3d.py
@@ -68,7 +68,8 @@ def get_edges_features_matrix_numba(elements, nodes):
                     ]
                 )
 
-    return edges_features_matrix, element_initial_volume
+    # TODO: #115 compute local_stifness_matrices which should be returned as 3rd value
+    return edges_features_matrix, element_initial_volume, None
 
 
 @numba.njit
diff --git a/conmech/dynamics/factory/dynamics_factory_method.py b/conmech/dynamics/factory/dynamics_factory_method.py
index da606485..3d4ffa6f 100644
--- a/conmech/dynamics/factory/dynamics_factory_method.py
+++ b/conmech/dynamics/factory/dynamics_factory_method.py
@@ -25,9 +25,11 @@ def get_basic_matrices(elements: np.ndarray, nodes: np.ndarray):
     dimension = len(elements[0]) - 1
     factory = get_factory(dimension)
 
-    edges_features_matrix, element_initial_volume = factory.get_edges_features_matrix(
-        elements, nodes
-    )
+    (
+        edges_features_matrix,
+        element_initial_volume,
+        local_stifness_matrices,
+    ) = factory.get_edges_features_matrix(elements, nodes)
 
     volume_at_nodes = edges_features_matrix[0]
     U = edges_features_matrix[1]
@@ -42,7 +44,7 @@ def get_basic_matrices(elements: np.ndarray, nodes: np.ndarray):
             for k in range(factory.dimension)
         ]
     )
-    return element_initial_volume, volume_at_nodes, U, V, W
+    return element_initial_volume, volume_at_nodes, U, V, W, local_stifness_matrices
 
 
 def get_dynamics(elements: np.ndarray, body_prop: BodyProperties, U, V, W):
diff --git a/conmech/simulations/problem_solver.py b/conmech/simulations/problem_solver.py
index 9e0f0096..81ee099b 100644
--- a/conmech/simulations/problem_solver.py
+++ b/conmech/simulations/problem_solver.py
@@ -109,6 +109,10 @@ def solving_method(self, value: str) -> None:
         self.__set_step_solver(value)
         self.__set_second_step_solver(value)
 
+    def refresh_solvers(self):
+        self.__set_step_solver(self.solving_method)
+        self.__set_second_step_solver(self.solving_method)
+
     def __set_step_solver(self, value):
         solver_class: Type[Solver] = SolversRegistry.get_by_name(
             solver_name=value, problem=self.problem
@@ -376,6 +380,12 @@ def solve(self, *, initial_displacement: Callable, verbose: bool = False, **kwar
         return state
 
 
+class NonHomogenousSolver(StaticSolver):
+    def update_density(self, density: np.ndarray):
+        self.body.reinitialize_matrices(density)
+        self.refresh_solvers()
+
+
 class QuasistaticRelaxation(ProblemSolver):
     def __init__(self, setup: RelaxationQuasistaticProblem, solving_method: str):
         """Solves general Contact Mechanics problem.
diff --git a/documentation/main.pdf b/documentation/main.pdf
new file mode 100644
index 00000000..08491fcf
Binary files /dev/null and b/documentation/main.pdf differ
diff --git a/documentation/main.tex b/documentation/main.tex
new file mode 100644
index 00000000..0fbc85ce
--- /dev/null
+++ b/documentation/main.tex
@@ -0,0 +1,29 @@
+\documentclass{article}
+
+\usepackage{url}
+\usepackage{amsmath}
+
+
+\title{Conmech documentation}
+
+\begin{document}
+
+\maketitle
+
+\subsection*{[LSM] Local stifness matrix}
+
+\[
+w_e^{kl} = \left[
+\int_{T_e} \frac{\partial \phi_i}{\partial x_k}\frac{\partial \phi_j}{\partial x_l}dx
+\right]_{i,j=1, \dots, N_e} \qquad \forall k,l \in \{1, 2\}
+\]
+
+\begin{tabular}{ll}
+Where & $e$ is the index of the element\\
+& $N_e$ is number of vertices of element $e$\\
+& $T_e$ is the surface of the mesh element $e$\\
+& $\phi_j$ is $j$-th base function of the element $e$
+\end{tabular}
+
+
+\end{document}
diff --git a/examples/example_density.py b/examples/example_density.py
new file mode 100644
index 00000000..d8f35e25
--- /dev/null
+++ b/examples/example_density.py
@@ -0,0 +1,63 @@
+"""
+Created at 02.09.2023
+"""
+from dataclasses import dataclass
+
+import numpy as np
+from conmech.helpers.config import Config
+from conmech.mesh.boundaries_description import BoundariesDescription
+from conmech.plotting.drawer import Drawer
+from conmech.scenarios.problems import StaticDisplacementProblem
+from conmech.simulations.problem_solver import NonHomogenousSolver
+
+from examples.p_slope_contact_law import make_slope_contact_law
+
+
+@dataclass
+class StaticSetup(StaticDisplacementProblem):
+    grid_height: ... = 1.0
+    elements_number: ... = (2, 5)
+    mu_coef: ... = 4
+    la_coef: ... = 4
+    contact_law: ... = make_slope_contact_law(slope=1)
+
+    @staticmethod
+    def inner_forces(x, t=None):
+        return np.array([0, 0])
+
+    @staticmethod
+    def outer_forces(x, t=None):
+        return np.array([0, 0.1]) if x[1] < 0.2 and x[0] > 2.2 else np.array([0, 0])
+
+    @staticmethod
+    def friction_bound(u_nu: float) -> float:
+        return 0
+
+    boundaries: ... = BoundariesDescription(
+        contact=lambda x: x[1] == 0 and x[0] < 0.2, dirichlet=lambda x: x[0] == 0
+    )
+
+
+def main(config: Config):
+    """
+    Entrypoint to example.
+
+    To see result of simulation you need to call from python `main(Config().init())`.
+    """
+    setup = StaticSetup(mesh_type="cross")
+    runner = NonHomogenousSolver(setup, "schur")
+
+    elem_centers = np.empty(shape=(len(runner.body.mesh.elements), 2))
+    for idx, elem in enumerate(runner.body.mesh.elements):
+        verts = runner.body.mesh.initial_nodes[elem]
+        elem_centers[idx] = np.sum(verts, axis=0) / len(elem)
+    density = np.asarray([1 if x[0] < 1 else 0.2 for x in elem_centers])
+
+    runner.update_density(density)
+
+    state = runner.solve(verbose=True, initial_displacement=setup.initial_displacement)
+    Drawer(state=state, config=config).draw(show=config.show, save=config.save)
+
+
+if __name__ == "__main__":
+    main(Config().init())
diff --git a/tests/test_conmech/unit_test/test_matrices.py b/tests/test_conmech/unit_test/test_matrices.py
index a78c5a84..64c3497a 100644
--- a/tests/test_conmech/unit_test/test_matrices.py
+++ b/tests/test_conmech/unit_test/test_matrices.py
@@ -22,7 +22,9 @@ def test_matrices_2d_integrals():
     )
 
     # Act
-    edges_features_matrix, element_initial_volume = sut_2d(elements=elements, nodes=initial_nodes)
+    edges_features_matrix, element_initial_volume, _ = sut_2d(
+        elements=elements, nodes=initial_nodes
+    )
 
     # Assert
     np.testing.assert_allclose(element_initial_volume.sum(), area)
@@ -50,7 +52,9 @@ def test_matrices_3d_integrals():
     )
 
     # Act
-    edges_features_matrix, element_initial_volume = sut_3d(elements=elements, nodes=initial_nodes)
+    edges_features_matrix, element_initial_volume, _ = sut_3d(
+        elements=elements, nodes=initial_nodes
+    )
 
     # Assert
     np.testing.assert_allclose(element_initial_volume.sum(), 1)