diff --git a/doc/chap_examples.rst b/doc/chap_examples.rst
index 2265c71b..c8b85f22 100644
--- a/doc/chap_examples.rst
+++ b/doc/chap_examples.rst
@@ -634,6 +634,26 @@ Defining :math:`x_1 = x, x_2=\dot{x}_1, x_3 = \theta` and :math:`x_4=\dot{x}_3`,
 
 Consider some moment of time :math:`t_1` and final position :math:`x_1(t_1) = x_1, x_2(t_1) = 0, x_3(t_1) = \frac{\pi}{2}, x_4(t_1) = 0`. It is required to calculate the backward reachability sets (tube) for the linearized system :eq:`invpendls` emanating from the given final position and project it onto :math:`(x_1, x_3)` subspace. It is also required to identify whether it’s possible to reach the final position from a given initial position :math:`x_1(t_1) = x_1, x_2(t_1) = v_1, x_3(t_1) = \theta_1, x_4(t_1) = \omega_1` using some admissible control function.
 
+.. literalinclude:: ../products/+elltool/+doc/+snip/s_chapter06_section05_snippet01.m
+   :language: matlab
+   :linenos:
+
+.. _brset:
+
+.. figure:: /pic/chapter06_section05_set1.png
+   :alt: set1
+   :width: 50 %
+
+   Backward reachability set
+
+.. _brtube:
+
+.. figure:: /pic/chapter06_section05_tube1.png
+   :alt: tube1
+   :width: 50 %
+
+   Backward reachability tube
+
 .. raw:: html
 	<h2>References</h2>
 
diff --git a/doc/pic/chapter06_section05_set1.png b/doc/pic/chapter06_section05_set1.png
new file mode 100644
index 00000000..e962abcf
Binary files /dev/null and b/doc/pic/chapter06_section05_set1.png differ
diff --git a/doc/pic/chapter06_section05_tube1.png b/doc/pic/chapter06_section05_tube1.png
new file mode 100644
index 00000000..e92a3bb3
Binary files /dev/null and b/doc/pic/chapter06_section05_tube1.png differ
diff --git a/products/+elltool/+doc/+snip/s_chapter06_section05_snippet01.m b/products/+elltool/+doc/+snip/s_chapter06_section05_snippet01.m
new file mode 100644
index 00000000..93ea6a10
--- /dev/null
+++ b/products/+elltool/+doc/+snip/s_chapter06_section05_snippet01.m
@@ -0,0 +1,48 @@
+alpha = 10;
+k1 = 50;
+k2 = 500;
+m1 = 150;
+m2 = 100;
+L = 5;
+g = 9.8;
+Jc = 100;
+
+x1 = 1;
+v1 = 0;
+theta1 = pi/2;
+omega1 = 0;
+t1 = 1;
+
+d = Jc + m2*L^2/4;
+% define matrices aMat, bMat, and control bounds uBoundsEllObj:
+aMat = [0, 1, 0, 0; 0, -k1/m1, 0, 0; 0, 0, 0, 1; 0, 0, m2*L*g/2/d, -k2*L/d];
+bMat = [0; 1/m1; 0; 0];
+uBoundsEllObj = alpha * ell_unitball(1);
+%define disturbance:
+gMat = [0; 0; 0; -m2*L*g*pi/4/d];
+vEllObj = ellipsoid(1, 0); %known disturbance
+
+%linear system
+lsys = elltool.linsys.LinSysContinuous(aMat, bMat, uBoundsEllObj,...
+    gMat, vEllObj);
+timeVec = [t1, 0];
+% initial directions (some random vectors in R^4):
+dirsMat = [1 0 1 0; 1 -1 0 0; 0 -1 0 1; 1 1 -1 1; -1 1 1 0; -2 0 1 1].';
+% x1EllObj = ellipsoid([x1; v1; theta1; omega1], eps*eye(4)); %known final point
+x1EllObj = ellipsoid([x1; v1; theta1; omega1], zeros(4)); %known final point
+
+
+%backward reach set
+brsObj = elltool.reach.ReachContinuous(lsys, x1EllObj, dirsMat, timeVec,...
+    'isRegEnabled', true, 'isJustCheck', false, 'regTol', 1e-4,...
+    'absTol', 1e-5, 'relTol', 1e-4);
+
+basisMat = [1 0 0 0; 0 0 1 0].'; % orthogonal basis of (x1, x3) subspace
+psObj = brsObj.projection(basisMat); % reach set projection
+
+% plot projection of reach set internal approximation:
+psObj.plotByIa('b');
+
+%plot backward reach set approximation at time t=0
+psCutObj = psObj.cut(0);
+psCutObj.plotByIa('r');