Classification and stability of fixed points of 2-D linear maps

CASE 1 a: Distinct real eigenvalues |λ|’s < 1 – stable NODE, (sink)

Clear[map,icoord,A]
imsize=420;xlim=2;ylim=2;rectangleE=1.7;
frame=15;
points=1200;(*points per frame*)
icoord=RandomReal[{-rectangleE,rectangleE},{points,2}];(*IC: random points in a square*)
A=(({{-0.7, 0}, {0, 0.9}}));(*system matrix*)
ev=Eigenvectors[A];
map[ic_,a_]:={(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[1,1]],(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[2,1]]}
sol=NestList[
Table[map[#,i],{i,1,points,1}]
&,icoord,frame];

ListAnimate[
Table[Show[
ListPlot[sol[[i]],PlotRange->{{-xlim,xlim},{-ylim,ylim}},Frame-> True,PlotStyle->Black,PlotMarkers->{"•",6},FrameLabel->{x,y},AspectRatio->0.9,LabelStyle->18,ImageSize->imsize],
Graphics[{Thick,Red,Arrow[{{0,0},#}]&/@ev}],
Graphics[{PointSize[0.03],Blue,Point[{0,0}]}]],
{i,1,frame}],
AnimationRunning->False,SaveDefinitions->True]

StringForm["Eigenvalues are ``, minus forces the iterates to jump around the corresponding manifold.",Eigenvalues[A]] 
StringForm["Determinant is ``, absolute value determines the value of the resulting area.",Det[A]] 
StringForm["Matrix trace is ``",Tr[A]] 
StringForm["Eigenvectors are ``",ev]

CASE 1 a: Distinct real eigenvalues |λ|’s > 1  – unstable NODE, (source)

Clear[map,icoord,A]
imsize=420;xlim=2;ylim=2;rectangleE=0.6;
frame=15;
points=1200;(*points per frame*)
icoord=RandomReal[{-rectangleE,rectangleE},{points,2}];(*IC: random points in a square*)
A=(({{0.0, 1}, {1.15, 0.01}}));(*system matrix*)
ev=Eigenvectors[A];

map[ic_,a_]:={(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[1,1]],(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[2,1]]}
sol=NestList[
Table[map[#,i],{i,1,points,1}]
&,icoord,frame];

ListAnimate[
Table[Show[
ListPlot[sol[[i]],PlotRange->{{-xlim,xlim},{-ylim,ylim}},Frame-> True,PlotStyle->Black,PlotMarkers->{"•",6},FrameLabel->{x,y},AspectRatio->0.9,LabelStyle->18,ImageSize->imsize],Graphics[{Thick,Red,Arrow[{{0,0},#}]&/@ev}],
Graphics[{PointSize[0.03],Blue,Point[{0,0}]}]],
{i,1,frame}],
AnimationRunning->False,SaveDefinitions->True]
StringForm["Eigenvalues are ``",Eigenvalues[A]] 
StringForm["Determinant is ``, absolute value determines the value of the resulting area.",Det[A]] 
StringForm["Matrix trace is ``",Tr[A]] 
StringForm["Eigenvectors are ``",ev]

CASE 1 b: Distinct real eigenvalues  |(Subscript[λ, 1])| > 1 > |(Subscript[λ, 2])| or |(Subscript[λ, 1])| < 1 < |(Subscript[λ, 2])| – SADDLE

One manifold is attracting and another is repelling

Clear[map,icoord,A]
imsize=420;xlim=2;ylim=2;rectangleE=1;
frame=15;
points=1200;(*points per frame*)
icoord=RandomReal[{-rectangleE,rectangleE},{points,2}];(*IC: random points in a square*)
A=(({{1.05, 0.4}, {0.1, -0.8}}));(*system matrix*)

ev=Eigenvectors[A];
map[ic_,a_]:={(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[1,1]],(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[2,1]]}
sol=NestList[
Table[map[#,i],{i,1,points,1}]
&,icoord,frame];

ListAnimate[
Table[Show[
ListPlot[sol[[i]],PlotRange->{{-xlim,xlim},{-ylim,ylim}},Frame-> True,PlotStyle->Black,PlotMarkers->{"•",6},FrameLabel->{x,y},AspectRatio->0.9,LabelStyle->18,ImageSize->imsize],Graphics[{Thick,Red,Arrow[{{0,0},#}]&/@ev}],
Graphics[{PointSize[0.03],Blue,Point[{0,0}]}]],
{i,1,frame}],
AnimationRunning->False,SaveDefinitions->True]
StringForm["Eigenvalues are `` (minus sign: jumps around corresponding manifold)",Eigenvalues[A]] 
StringForm["Determinant is ``, absolute value determines the value of the resulting area.",Det[A]] 
StringForm["Matrix trace is ``",Tr[A]] 
StringForm["Eigenvectors are ``",ev]

CASE 1 c: Distinct real eigenvalues  |(Subscript[λ, 1])| > 1 > |(Subscript[λ, 2])| or |(Subscript[λ, 1])| < 1 < |(Subscript[λ, 2])| – SADDLE

One manifold is attracting and another is repelling

Clear[map,icoord,A]
imsize=420;xlim=2;ylim=2;rectangleE=1;
frame=15;
points=1200;(*points per frame*)
icoord=RandomReal[{-rectangleE,rectangleE},{points,2}];(*IC: random points in a square*)
A=(({{-1.05, 0.4}, {0.1, -0.8}}));(*system matrix*)

ev=Eigenvectors[A];
map[ic_,a_]:={(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[1,1]],(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[2,1]]}
sol=NestList[
Table[map[#,i],{i,1,points,1}]
&,icoord,frame];

ListAnimate[
Table[Show[
ListPlot[sol[[i]],PlotRange->{{-xlim,xlim},{-ylim,ylim}},Frame-> True,PlotStyle->Black,PlotMarkers->{"•",6},FrameLabel->{x,y},AspectRatio->0.9,LabelStyle->18,ImageSize->imsize],Graphics[{Thick,Red,Arrow[{{0,0},#}]&/@ev}],
Graphics[{PointSize[0.03],Blue,Point[{0,0}]}]],
{i,1,frame}],
AnimationRunning->False,SaveDefinitions->True]
StringForm["Eigenvalues are `` (minus sign: jumps around corresponding manifold)",Eigenvalues[A]] 
StringForm["Determinant is ``, absolute value determines the value of the resulting area.",Det[A]] 
StringForm["Matrix trace is ``",Tr[A]] 
StringForm["Eigenvectors are ``",ev]

CASE 2 a: Repeated real eigenvalues  (Subscript[λ, 1])=(Subscript[λ, 2]) – DEGENERATE NODE, stable if |λ| <1 and Δ<1 (sink) and unstable if |λ| >1 and Δ>1 (source)

Only one eigen direction

Clear[map,icoord,A]
imsize=420;xlim=2;ylim=2;rectangleE=0.5;
frame=15;
points=1200;(*points per frame*)
icoord=RandomReal[{-rectangleE,rectangleE},{points,2}];(*IC: random points in a square*)
A=(({{1.1, 1}, {0, 1.1}}));(*system matrix (({{-1, 2}, {0, -1}}))*)

ev=Eigenvectors[A];
map[ic_,a_]:={(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[1,1]],(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[2,1]]}
sol=NestList[
Table[map[#,i],{i,1,points,1}]
&,icoord,frame];

ListAnimate[
Table[Show[
ListPlot[sol[[i]],PlotRange->{{-xlim,xlim},{-ylim,ylim}},Frame-> True,PlotStyle->Black,PlotMarkers->{"•",6},FrameLabel->{x,y},AspectRatio->0.9,LabelStyle->18,ImageSize->imsize],Graphics[{Thick,Red,Arrow[{{0,0},#}]&/@ev}],
Graphics[{PointSize[0.03],Blue,Point[{0,0}]}]],
{i,1,frame}],
AnimationRunning->False,SaveDefinitions->True]
StringForm["Eigenvalues are ``",Eigenvalues[A]] 
StringForm["Determinant is ``, absolute value determines the value of the resulting area.",Det[A]] 
StringForm["Matrix trace is ``",Tr[A]] 
StringForm["Eigenvectors are ``",ev] 
StringForm["(τ^2)-4Δ = `` (sits on the parabola)",Tr[A]^2-4 Det[A]]

CASE 2 a: Repeated real eigenvalues  (Subscript[λ, 1])=(Subscript[λ, 2]) – DEGENERATE NODE, stable if |λ| <1 and Δ<1 (sink) and unstable if |λ| >1 and Δ>1 (source)

Only one eigen direction

Clear[map,icoord,A]
imsize=420;xlim=2;ylim=2;rectangleE=1.5;
frame=15;
points=1200;(*points per frame*)
icoord=RandomReal[{-rectangleE,rectangleE},{points,2}];(*IC: random points in a square*)
A=(({{0.7, 1}, {0, 0.7}}));(*system matrix (({{-1, 2}, {0, -1}}))*)

ev=Eigenvectors[A];
map[ic_,a_]:={(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[1,1]],(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[2,1]]}
sol=NestList[
Table[map[#,i],{i,1,points,1}]
&,icoord,frame];

ListAnimate[
Table[Show[
ListPlot[sol[[i]],PlotRange->{{-xlim,xlim},{-ylim,ylim}},Frame-> True,PlotStyle->Black,PlotMarkers->{"•",6},FrameLabel->{x,y},AspectRatio->0.9,LabelStyle->18,ImageSize->imsize],Graphics[{Thick,Red,Arrow[{{0,0},#}]&/@ev}],
Graphics[{PointSize[0.03],Blue,Point[{0,0}]}]],
{i,1,frame}],
AnimationRunning->False,SaveDefinitions->True]
StringForm["Eigenvalues are ``",Eigenvalues[A]] 
StringForm["Determinant is ``, absolute value determines the value of the resulting area.",Det[A]] 
StringForm["Matrix trace is ``",Tr[A]] 
StringForm["Eigenvectors are ``",ev] 
StringForm["(τ^2)-4Δ = `` (sits on the parabola)",Tr[A]^2-4 Det[A]]

CASE 2 b: Repeated real eigenvalues  (Subscript[λ, 1])=(Subscript[λ, 2])  – degenerate node - STAR, stable if |λ|’s <1 and Δ<1 (sink) and unstable if |λ|’s >1 and Δ>1 (source)

All directions are eigen directions

Clear[map,icoord,A]
imsize=420;xlim=2;ylim=2;rectangleE=1.5;
frame=15;
points=1200;(*points per frame*)
icoord=RandomReal[{-rectangleE,rectangleE},{points,2}];(*IC: random points in a square*)
A=(({{0.8, 0}, {0, 0.8}}));(*system matrix*)

ev=Eigenvectors[A];
map[ic_,a_]:={(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[1,1]],(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[2,1]]}
sol=NestList[
Table[map[#,i],{i,1,points,1}]
&,icoord,frame];

ListAnimate[
Table[Show[
ListPlot[sol[[i]],PlotRange->{{-xlim,xlim},{-ylim,ylim}},Frame-> True,PlotStyle->Black,PlotMarkers->{"•",6},FrameLabel->{x,y},AspectRatio->0.9,LabelStyle->18,ImageSize->imsize],Graphics[{Thick,Red,Arrow[{{0,0},#}]&/@ev}],
Graphics[{PointSize[0.03],Blue,Point[{0,0}]}]],
{i,1,frame}],
AnimationRunning->False,SaveDefinitions->True]
StringForm["Eigenvalues are ``, if eigenvaluea are negative points jump over the origin",Eigenvalues[A]] 
StringForm["Determinant is ``, absolute value determines the value of the resulting area.",Det[A]] 
StringForm["Matrix trace is ``",Tr[A]] 
StringForm["Eigenvectors are ``",ev] 
StringForm["(τ^2)-4Δ = `` (sits on the parabola, see the classification chart)",Tr[A]^2-4 Det[A]]

CASE 2 b: Repeated real eigenvalues  (Subscript[λ, 1])=(Subscript[λ, 2])  – degenerate node - STAR, stable if |λ|’s <1 and Δ<1 (sink) and unstable if |λ|’s >1 and Δ>1 (source)

All directions are eigen directions

Clear[map,icoord,A]
imsize=420;xlim=2;ylim=2;rectangleE=0.5;
frame=15;
points=1200;(*points per frame*)
icoord=RandomReal[{-rectangleE,rectangleE},{points,2}];(*IC: random points in a square*)
A=(({{-1.1, 0}, {0, -1.1}}));(*system matrix*)

ev=Eigenvectors[A];
map[ic_,a_]:={(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[1,1]],(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[2,1]]}
sol=NestList[
Table[map[#,i],{i,1,points,1}]
&,icoord,frame];

ListAnimate[
Table[Show[
ListPlot[sol[[i]],PlotRange->{{-xlim,xlim},{-ylim,ylim}},Frame-> True,PlotStyle->Black,PlotMarkers->{"•",6},FrameLabel->{x,y},AspectRatio->0.9,LabelStyle->18,ImageSize->imsize],Graphics[{Thick,Red,Arrow[{{0,0},#}]&/@ev}],
Graphics[{PointSize[0.03],Blue,Point[{0,0}]}]],
{i,1,frame}],
AnimationRunning->False,SaveDefinitions->True]
StringForm["Eigenvalues are ``, if eigenvaluea are negative points jump over the origin",Eigenvalues[A]] 
StringForm["Determinant is ``, absolute value determines the value of the resulting area.",Det[A]] 
StringForm["Matrix trace is ``",Tr[A]] 
StringForm["Eigenvectors are ``",ev] 
StringForm["(τ^2)-4Δ = `` (sits on the parabola, see the classification chart)",Tr[A]^2-4 Det[A]]

CASE 3 a: Complex eigenvalues – Linear rotation, stable if (Sqrt[((Re[(Subscript[λ, i])])^2)+((Im[(Subscript[λ, i])])^2)])<1  and Δ<1 (sink) and unstable if (Sqrt[((Re[(Subscript[λ, i])])^2)+((Im[(Subscript[λ, i])])^2)])>1  and Δ>1 (source)

Clear[map,icoord,A]
imsize=420;xlim=2;ylim=2;rectangleE=1.5;
frame=15;
points=1200;(*points per frame*)
icoord=RandomReal[{-rectangleE,rectangleE},{points,2}];(*IC: random points in a square*)

A=(({{0.1, -0.8}, {0.8, 0.1}}));(*system matrix*)

ev=Eigenvectors[A];
map[ic_,a_]:={(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[1,1]],(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[2,1]]}
sol=NestList[
Table[map[#,i],{i,1,points,1}]
&,icoord,frame];

ListAnimate[
Table[Show[
ListPlot[sol[[i]],PlotRange->{{-xlim,xlim},{-ylim,ylim}},Frame-> True,PlotStyle->Black,PlotMarkers->{"•",6},FrameLabel->{x,y},AspectRatio->0.9,LabelStyle->18,ImageSize->imsize],Graphics[{PointSize[0.03],Blue,Point[{0,0}]}]],
{i,1,frame}],
AnimationRunning->False,SaveDefinitions->True]
StringForm["Eigenvalues are ``",Eigenvalues[A]] 
StringForm["Determinant is ``, absolute value determines the value of the resulting area.",Det[A]] 
StringForm["Matrix trace is ``",Tr[A]] 
StringForm["Eigenvectors are ``",ev] 
StringForm["(τ^2)-4Δ = `` (sits under the parabola, (τ^2)-4Δ < 0)",Tr[A]^2-4 Det[A]] 
StringForm["(Sqrt[((Re[(Subscript[λ, 1])])^2) + ((Im[(Subscript[λ, 1])])^2)])= `` (stability condition)",N[(Sqrt[((Re[Eigenvalues[A][[1]]])^2)+((Im[Eigenvalues[A][[1]]])^2)])]] 
StringForm["(Sqrt[((Re[(Subscript[λ, 2])])^2) + ((Im[(Subscript[λ, 2])])^2)])= `` (stability condition)",N[(Sqrt[((Re[Eigenvalues[A][[2]]])^2)+((Im[Eigenvalues[A][[2]]])^2)])]]

CASE 3 a: Complex eigenvalues – Linear rotation, stable if (Sqrt[((Re[(Subscript[λ, i])])^2)+((Im[(Subscript[λ, i])])^2)])<1  and Δ<1 (sink) and unstable if (Sqrt[((Re[(Subscript[λ, i])])^2)+((Im[(Subscript[λ, i])])^2)])>1  and Δ>1 (source)

Clear[map,icoord,A]
imsize=420;xlim=2;ylim=2;rectangleE=0.5;
frame=15;
points=1200;(*points per frame*)
icoord=RandomReal[{-rectangleE,rectangleE},{points,2}];(*IC: random points in a square*)

A=(({{0.2, -1.1}, {1.1, 0.2}}));(*system matrix*)

ev=Eigenvectors[A];
map[ic_,a_]:={(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[1,1]],(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[2,1]]}
sol=NestList[
Table[map[#,i],{i,1,points,1}]
&,icoord,frame];

ListAnimate[
Table[Show[
ListPlot[sol[[i]],PlotRange->{{-xlim,xlim},{-ylim,ylim}},Frame-> True,PlotStyle->Black,PlotMarkers->{"•",6},FrameLabel->{x,y},AspectRatio->0.9,LabelStyle->18,ImageSize->imsize],Graphics[{PointSize[0.03],Blue,Point[{0,0}]}]],
{i,1,frame}],
AnimationRunning->False,SaveDefinitions->True]
StringForm["Eigenvalues are ``",Eigenvalues[A]] 
StringForm["Determinant is ``, absolute value determines the value of the resulting area.",Det[A]] 
StringForm["Matrix trace is ``",Tr[A]] 
StringForm["Eigenvectors are ``",ev] 
StringForm["(τ^2)-4Δ = `` (sits under the parabola, (τ^2)-4Δ < 0)",Tr[A]^2-4 Det[A]] 
StringForm["(Sqrt[((Re[(Subscript[λ, 1])])^2) + ((Im[(Subscript[λ, 1])])^2)])= `` (stability condition)",N[(Sqrt[((Re[Eigenvalues[A][[1]]])^2)+((Im[Eigenvalues[A][[1]]])^2)])]] 
StringForm["(Sqrt[((Re[(Subscript[λ, 2])])^2) + ((Im[(Subscript[λ, 2])])^2)])= `` (stability condition)",N[(Sqrt[((Re[Eigenvalues[A][[2]]])^2)+((Im[Eigenvalues[A][[2]]])^2)])]]

CASE 3 a: Complex eigenvalues – Linear rotation, stable if (Sqrt[((Re[(Subscript[λ, i])])^2)+((Im[(Subscript[λ, i])])^2)])<1 and Δ<1 (sink) and unstable if (Sqrt[((Re[(Subscript[λ, i])])^2)+((Im[(Subscript[λ, i])])^2)])>1  and Δ>1 (source)

Clear[map,icoord,A]
imsize=420;xlim=2;ylim=2;rectangleE=0.5;
frame=15;
points=1200;(*points per frame*)
icoord=RandomReal[{-rectangleE,rectangleE},{points,2}];(*IC: random points in a square*)

A=(({{0.7, 0.8}, {-0.8, 0.7}}));(*system matrix*)

ev=Eigenvectors[A];
map[ic_,a_]:={(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[1,1]],(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[2,1]]}
sol=NestList[
Table[map[#,i],{i,1,points,1}]
&,icoord,frame];

ListAnimate[
Table[Show[
ListPlot[sol[[i]],PlotRange->{{-xlim,xlim},{-ylim,ylim}},Frame-> True,PlotStyle->Black,PlotMarkers->{"•",6},FrameLabel->{x,y},AspectRatio->0.9,LabelStyle->18,ImageSize->imsize],Graphics[{PointSize[0.03],Blue,Point[{0,0}]}]],
{i,1,frame}],
AnimationRunning->False,SaveDefinitions->True]
StringForm["Eigenvalues are ``",Eigenvalues[A]] 
StringForm["Determinant is ``, absolute value determines the value of the resulting area.",Det[A]] 
StringForm["Matrix trace is ``",Tr[A]] 
StringForm["Eigenvectors are ``",ev] 
StringForm["(τ^2)-4Δ = `` (sits under the parabola, (τ^2)-4Δ < 0)",Tr[A]^2-4 Det[A]] 
StringForm["(Sqrt[((Re[(Subscript[λ, 1])])^2) + ((Im[(Subscript[λ, 1])])^2)])= `` (stability condition)",N[(Sqrt[((Re[Eigenvalues[A][[1]]])^2)+((Im[Eigenvalues[A][[1]]])^2)])]] 
StringForm["(Sqrt[((Re[(Subscript[λ, 2])])^2) + ((Im[(Subscript[λ, 2])])^2)])= `` (stability condition)",N[(Sqrt[((Re[Eigenvalues[A][[2]]])^2)+((Im[Eigenvalues[A][[2]]])^2)])]]

CASE 3 b: Pure imaginary complex eigenvalues – Linear rotation where (Sqrt[((Re[(Subscript[λ, i])])^2)+((Im[(Subscript[λ, i])])^2)])=1 and Δ=1

The area remains unchanged, Δ = det(A) = 1

Clear[map,icoord,A]
imsize=420;xlim=2;ylim=2;rectangleE=1.4;
frame=15;
points=1000;(*points per frame*)
icoord=RandomReal[{-rectangleE,rectangleE},{points,2}];(*IC: random points in a square*)

θ=(π/12); (*rotation angle try general one e.g. π/11*)
A=(({{Cos[θ], -Sin[θ]}, {Sin[θ], Cos[θ]}}));(*system matrix, general incudes case III, a)*)

ev=Eigenvectors[A];
map[ic_,a_]:={(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[1,1]],(A.(({{ic[[a,1]]}, {ic[[a,2]]}})))[[2,1]]}
sol=NestList[
Table[map[#,i],{i,1,points,1}]
&,icoord,frame];

ListAnimate[
Table[Show[
ListPlot[sol[[i]],PlotRange->{{-xlim,xlim},{-ylim,ylim}},Frame-> True,PlotStyle->Black,PlotMarkers->{"•",6},FrameLabel->{x,y},AspectRatio->0.9,LabelStyle->18,ImageSize->imsize],Graphics[{PointSize[0.03],Blue,Point[{0,0}]}]],
{i,1,frame}],
AnimationRunning->False,SaveDefinitions->True]
StringForm["Eigenvalues are ``",N[Eigenvalues[A]]] 
StringForm["Determinant is ``, absolute value determines the value of the resulting area.",Det[A]] 
StringForm["Matrix trace is ``",Tr[A]] 
StringForm["Eigenvectors are ``",ev] 
StringForm["(τ^2)-4Δ = `` (sits under the parabola, (τ^2)-4Δ < 0)",Tr[A]^2-4 Det[A]] 
StringForm["(Sqrt[((Re[(Subscript[λ, 1])])^2) + ((Im[(Subscript[λ, 1])])^2)])= `` (stability condition)",N[(Sqrt[((Re[Eigenvalues[A][[1]]])^2)+((Im[Eigenvalues[A][[1]]])^2)])]] 
StringForm["(Sqrt[((Re[(Subscript[λ, 2])])^2) + ((Im[(Subscript[λ, 2])])^2)])= `` (stability condition)",N[(Sqrt[((Re[Eigenvalues[A][[2]]])^2)+((Im[Eigenvalues[A][[2]]])^2)])]]

Creation of this notebook was supported by IT Academy program of Information Technology Foundation for Education (HITSA).