diff --git a/.gitignore b/.gitignore
index 961f132..c26000a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,5 @@
+# LaTeX
+build
+
 # MATLAB
 *.asv
diff --git a/Direct Methods for Solving Linear Systems/Alg1.m b/01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg1.m
similarity index 100%
rename from Direct Methods for Solving Linear Systems/Alg1.m
rename to 01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg1.m
diff --git a/Direct Methods for Solving Linear Systems/Alg10.m b/01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg10.m
similarity index 100%
rename from Direct Methods for Solving Linear Systems/Alg10.m
rename to 01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg10.m
diff --git a/Direct Methods for Solving Linear Systems/Alg11.m b/01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg11.m
similarity index 100%
rename from Direct Methods for Solving Linear Systems/Alg11.m
rename to 01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg11.m
diff --git a/Direct Methods for Solving Linear Systems/Alg13.m b/01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg13.m
similarity index 100%
rename from Direct Methods for Solving Linear Systems/Alg13.m
rename to 01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg13.m
diff --git a/Direct Methods for Solving Linear Systems/Alg2.m b/01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg2.m
similarity index 100%
rename from Direct Methods for Solving Linear Systems/Alg2.m
rename to 01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg2.m
diff --git a/Direct Methods for Solving Linear Systems/Alg3.m b/01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg3.m
similarity index 100%
rename from Direct Methods for Solving Linear Systems/Alg3.m
rename to 01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg3.m
diff --git a/Direct Methods for Solving Linear Systems/Alg4.m b/01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg4.m
similarity index 100%
rename from Direct Methods for Solving Linear Systems/Alg4.m
rename to 01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg4.m
diff --git a/Direct Methods for Solving Linear Systems/Alg5.m b/01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg5.m
similarity index 100%
rename from Direct Methods for Solving Linear Systems/Alg5.m
rename to 01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg5.m
diff --git a/Direct Methods for Solving Linear Systems/Alg6_RREF.m b/01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg6_RREF.m
similarity index 100%
rename from Direct Methods for Solving Linear Systems/Alg6_RREF.m
rename to 01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg6_RREF.m
diff --git a/Direct Methods for Solving Linear Systems/Alg7.m b/01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg7.m
similarity index 100%
rename from Direct Methods for Solving Linear Systems/Alg7.m
rename to 01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg7.m
diff --git a/Direct Methods for Solving Linear Systems/Alg8.m b/01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg8.m
similarity index 100%
rename from Direct Methods for Solving Linear Systems/Alg8.m
rename to 01_Direct-Methods-for-Solving-Linear-Systems/Code/Alg8.m
diff --git a/Direct Methods for Solving Linear Systems/householder.m b/01_Direct-Methods-for-Solving-Linear-Systems/Code/householder.m
similarity index 100%
rename from Direct Methods for Solving Linear Systems/householder.m
rename to 01_Direct-Methods-for-Solving-Linear-Systems/Code/householder.m
diff --git a/Direct Methods for Solving Linear Systems/main.m b/01_Direct-Methods-for-Solving-Linear-Systems/Code/main.m
similarity index 99%
rename from Direct Methods for Solving Linear Systems/main.m
rename to 01_Direct-Methods-for-Solving-Linear-Systems/Code/main.m
index 502ca22..3afcbed 100644
--- a/Direct Methods for Solving Linear Systems/main.m	
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Code/main.m
@@ -23,8 +23,6 @@ b = P*b;
 y = Alg3(L, b);     % Forward substitution
 x = Q*Alg4(U, y)	% Backward substitution
 
-
-
 %% Problem 4
 
 A = [0.835, 0.667;
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/README.md b/01_Direct-Methods-for-Solving-Linear-Systems/README.md
new file mode 100644
index 0000000..23f7b7f
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/README.md
@@ -0,0 +1,10 @@
+## Lab 1 - Direct Methods for Solving Linear Systems
+
+
+
+### Bibliography
+
+References in the code comments are represented by `[number]`:
+
+  - 1
+  - 2 – Golub, van Loan, Matrix Computations, 3rd edition
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/IMG_0216.jpg b/01_Direct-Methods-for-Solving-Linear-Systems/Report/IMG_0216.jpg
new file mode 100644
index 0000000..a4ec0a7
Binary files /dev/null and b/01_Direct-Methods-for-Solving-Linear-Systems/Report/IMG_0216.jpg differ
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg1.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg1.m
new file mode 100644
index 0000000..2c66149
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg1.m
@@ -0,0 +1,26 @@
+function A = Alg1(A)
+% Algorithm 1: Outer Product Gaussian Elimination
+% Performs a gaussian eliminaion on a square matrix A.
+
+[m, n] = size(A);
+
+if m ~= n
+    error('Matrix is not square!')
+end
+
+for k = 1:n-1
+    if det(A(1:k, 1:k)) < eps
+        error('Matrix is not nonsingular!')
+    end
+end
+
+% The following algorithm is based on the Algrotihm 3.2.1 from [2].
+
+for k = 1 : m-1
+    rows = k + 1 : m;
+    A(rows, k) = A(rows, k)/A(k, k);
+    A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows);
+end
+ 
+end
+
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg10.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg10.m
new file mode 100644
index 0000000..e17c80e
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg10.m
@@ -0,0 +1,24 @@
+function A = Alg10(A)
+% Algorithm 10: Cholesky (Banachiewicz) factorization.
+% Matrix A has to be symmetric and positive-definite:
+% all d(i, i) > 0 for A = LDL^T.
+% We are using Outer product Version (Golub, Load, Alg. 4.2.2), which
+% computes lower triangular G, such that A = G*G^T
+
+[m, n] = size(A);
+
+if m ~= n
+    error('Matrix is not square!')
+end
+
+for k = 1:m
+   A(k,k) = sqrt(A(k,k));
+   A(k+1:n, k) = A(k+1:n, k)/A(k,k);
+   for j = k+1:n
+       A(j:n, j) = A(j:n, j) - A(j:n, k)*A(j, k);
+   end
+end
+
+A = tril(A);
+
+end
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg11.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg11.m
new file mode 100644
index 0000000..6c09fda
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg11.m
@@ -0,0 +1,22 @@
+function [Q R] = Alg11(A)
+% Algorithm 11: QR factorization via Householder transformation.
+
+[m, n] = size(A);
+
+if ~ (m>=n)
+   error("m has to be greater or equal n!") 
+end
+
+for j = 1:n
+    [v, beta] = householder(A(j:m, j))
+    A(j:m, j:n) = (eye(m-j+1)-beta*(v*v.'))*A(j:m, j:n)
+    if j < m
+        A(j+1:m, j) = v(2:m-j+1)
+    end
+end
+
+R = triu(A(1:n, 1:n))
+
+Q = tril(A, -1)
+
+end
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg12.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg12.m
new file mode 100644
index 0000000..0538680
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg12.m
@@ -0,0 +1,21 @@
+function [Q, R] = Alg12(A)
+% Algorithm 11 - The  QR  algorithm  by  the  Householder  transformation
+% [Q, R] = Alg12(A)
+% Based on [GVL96]: Algorithm 5.2.4
+[m, n] = size(A);
+
+if ~ (m >= n)
+   error("m has to be greater or equal n!") 
+end
+
+Q = eye(m, m);
+for j = 1:n
+    for i = m:-1:j+1
+        [c, s] = Givens(A(i-1, j), A(i, j));
+        A(i-1:i, j:n) = [c s; -s c]'*A(i-1:i, j:n);
+        Q(:, i-1:i) = Q(:, i-1:i)*[c s ; -s c];
+    end
+end
+R = A;
+  
+end
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg13.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg13.m
new file mode 100644
index 0000000..89d40be
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg13.m
@@ -0,0 +1,19 @@
+function [Q, R] = Alg13(A)
+% Algorithm 13: The QR algorithm by the Gram-Schmidt ortogonalization
+% [Q, R] = Alg13(A)
+% Based on [GVL96]: Section 5.2.8
+
+[m,n] = size(A);
+
+if ~ (m >= n)
+   error("m has to be greater or equal n!") 
+end
+
+R = zeros(n, n);
+Q = zeros(m, n);
+for k = 1:n
+    R(k, k) = norm(A(:, k));
+    Q(:, k) = A(:, k)/R(k, k);
+    R(k, k+1:n) = Q(:, k)'*A(:, k+1:n);
+    A(:, k+1:n) = A(:, k+1:n) - Q(:, k)*R(k, k+1:n);
+end
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg2.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg2.m
new file mode 100644
index 0000000..30be69d
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg2.m
@@ -0,0 +1,47 @@
+function [P, Q, L, U] = Alg2(A)
+% Algorithm 2: Gaussian Elimination with Complete Pivoting.
+% [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A)
+% computes the complete pivoting factorization PAQ = LU.
+
+[m, n] = size(A);
+
+if m ~= n
+    error('Matrix is not square!')
+end
+
+% p and q are permutation vectors - respectively rows and columns
+p = 1:m;
+q = 1:m;
+
+% The following algorithm is based on the Algrotihm 3.4.2 from [2].
+
+for k = 1 : m-1
+    i = k:m;
+    j = k:m;
+    % Find the maximum entry to be the next pivot
+    [max_val, rows_of_max_in_col] = max(abs(A(i, j)));
+    [~, max_col] = max(max_val);
+    max_row = rows_of_max_in_col(max_col);
+    % Assign value of mu and lambda in respect to the main matrix A
+    [mi, lm] = deal(max_row+k-1, max_col+k-1);
+    % Interchange the rows and columns of matrix A...
+    A([k mi], 1:m) = deal(A([mi k], 1:m));
+    A(1:m, [k lm]) = deal(A(1:m, [lm k]));
+    % ...and respective permutation vectors entries.
+    p([k, mi]) = p([mi, k]);
+    q([k, lm]) = q([lm, k]);
+    % Perform Gaussian elimination with the greatest pivot
+    if A(k, k) ~= 0
+        rows = k+1 : m;
+        A(rows, k) = A(rows, k)/A(k, k);
+        A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows);
+    end
+end
+
+I = eye(m);
+U = triu(A);
+L = tril(A, -1) + I;
+P = I(p, :);
+Q = I(:, q);
+
+end
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg3.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg3.m
new file mode 100644
index 0000000..17af4b5
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg3.m
@@ -0,0 +1,29 @@
+function b = Al3(L, b)
+% Algorithm 3: Forward Substitution
+% b = Alg3_forward_substitution(L, b) overwrites b with the solution to
+% Lx = b.
+
+[m, n] = size(L);
+
+if m ~= n
+    error('Matrix is not square!')
+end
+
+if length(b) ~= m
+    error('Vector b has wrong length!')
+end
+
+if det(L) < eps
+    error("Matrix is not nonsingular!")
+end
+
+% The following algorithm is based on the Algrotihm 3.1.1 from [2].
+
+% b(m, :) so that matrices are also accepted
+b(1, :) = b(1, :)/L(1,1);
+for i = 2:m
+    b(i, :) = (b(i, :) - L(i, 1:i-1)*b(1:i-1, :))/L(i, i);
+end
+
+end
+
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg4.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg4.m
new file mode 100644
index 0000000..5a6d538
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg4.m
@@ -0,0 +1,33 @@
+function b = Alg4(U,b)
+% Argorithm 4: Back Substitution 
+% b = Alg4_back_substitution(U,b) returns vetor b with solution to the
+% Ux = b.
+
+[m, n] = size(U);
+
+if U ~= triu(U)
+    error('Matrix U is not upper triangular!')
+end
+
+if m ~= n
+    error('Matrix is not squared!')
+end
+
+if length(b) ~= m
+    error('Vector b has wrong length!')
+end
+
+% if det(U) < eps
+%     error('Matrix is not nonsingular!')
+% end
+
+
+% The following algorithm is based on the Algrotihm 3.1.2 from [2].
+
+% b(m, :) so that matrices are also accepted
+b(m, :) = b(m, :)/U(m, m);
+for i = m-1:-1:1
+    b(i, :) = (b(i, :) - U(i, i+1 : m)*b(i+1 : m, :))/U(i, i);
+end
+
+end
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg5.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg5.m
new file mode 100644
index 0000000..6c83199
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg5.m
@@ -0,0 +1,21 @@
+function A = Alg5(A)
+% Algorithm 5: Gauss-Jordan Elimination
+% A = Alg5_gauss_jordan_elimination(A) performs Gauss-Jordan elimination
+% on an augmented matrix A.
+
+[m, ~] = size(A);
+
+for k = 1 : m
+    
+    row = A(k, :);
+    row = row/row(k);
+    A(k, :) = row;
+    
+    for i = 1 : m
+        if i ~= k && A(i, k) ~= 0
+            A(i, :) = A(i, :)-(A(i, k))*row;
+        end
+    end
+end
+
+end
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg6_RREF.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg6_RREF.m
new file mode 100644
index 0000000..c57bb83
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg6_RREF.m
@@ -0,0 +1,46 @@
+function A = Alg6_RREF(A)
+% Algorithm 6: Reduced Row Echelon Form (RREF)
+% A = Alg6_RREF(A) returns RREF of matrix A. 
+
+[m, n] = size(A);
+
+j = 0;
+
+for k = 1 : m
+    j = j + 1;
+    if j > n
+        break
+    end
+    % We want the left-most coefficient to be 1 (pivot)
+    row = A(k, :);
+    if row(k) == 0
+        j = j + 1;
+    end
+    row = row/row(j);
+    A(k, :) = row;
+    
+    for i = 1 : m
+        if i ~= k
+            A(i, :) = A(i, :)-(A(i, j))*row;
+        end
+    end
+    
+    for i = k + 1 : m
+       A(i:end, k+1:end); % Partial matrix (in which we are looking for non-zero pivots)
+       A(i:end, k+1); % Left-most column
+       if ~any(A(i:end, k+1)) % If the left-most column has only zeros check the next one
+           k = k + 1;
+       end
+       A(i:end, k+1:end);
+       if A(i, k+1) == 0
+           non_zero_row = find(A(i:end,k+1), 1);
+           if isempty(non_zero_row)
+               continue
+           end
+           A([i, i+non_zero_row-1], :) = deal(A([i+non_zero_row-1, i], :));
+       end
+    end
+end
+
+end
+
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg7.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg7.m
new file mode 100644
index 0000000..e3a214e
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg7.m
@@ -0,0 +1,35 @@
+function [L, U] = Alg7(A)
+% LU Factorization without pivoting
+% [L, U] = Alg7(A) decomposes matrix A using Crout's algorithm into
+% U - upper triangular matrix and L - unit lower triangular matrix
+% such that A = LU.
+
+[m, n] = size(A);
+
+% This solution is mathematically correct, however computationaly
+% inefficient.
+%
+% A = Alg1_outer_product_gaussian_elimination(A);
+% 
+% U = triu(A);
+% L = tril(A, -1) + eye(m);
+
+% Instead, we should use the Crout's algorithm
+
+L = zeros(m, n);
+U = eye(m, n);
+
+for k = 1:n
+    
+    for j = k : n
+        U(k, j) = A(k, j) - dot(L(k, 1:k-1), U(1:k-1, j));
+    end
+    
+    for i = k:n
+        L(i, k) = (A(i, k) - dot(L(i,1:k-1), U(1:k-1, k))) / U(k,k);
+    end
+    
+    
+end
+
+end
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg8.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg8.m
new file mode 100644
index 0000000..6dfb4d4
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg8.m
@@ -0,0 +1,34 @@
+function [L, U, P] = Alg8(A)
+% LU Factorization with partial pivoting
+% [L, U, P] = Alg8(A) decomposes matrix A using Crout's algorithm into
+% U - upper triangular matrix and L - unit lower triangular matrix
+% using partial pivoting, interchanging A's rows, such that AP = LU.
+
+[m, n] = size(A);
+
+L = zeros(m, n);
+U = eye(m, n);
+P = eye(m, n);
+
+for k = 1:m
+    if A(k, k) == 0
+        non_zero_col = find(A(k, :), 1);
+        A(:, [k non_zero_col]) = deal(A(:, [non_zero_col k]));
+        P(:, [k non_zero_col]) = deal(P(:, [non_zero_col k]));
+    end
+end
+
+for k = 1:n
+    
+    for j = k : n
+        U(k, j) = A(k, j) - dot(L(k, 1:k-1), U(1:k-1, j));
+    end
+    
+    for i = k:n
+        L(i, k) = (A(i, k) - dot(L(i,1:k-1), U(1:k-1, k))) / U(k,k);
+    end
+    
+    
+end
+
+end
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/bibliography.bib b/01_Direct-Methods-for-Solving-Linear-Systems/Report/bibliography.bib
new file mode 100644
index 0000000..2aceb78
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/bibliography.bib
@@ -0,0 +1,13 @@
+@misc{Zdunek,
+    author = {Zdunek, Rafał},
+    title = {Lecture notes in {Numerical Algorithms}},
+    year = {2021}
+}
+
+@Book{GoluVanl96,
+  Title                    = {Matrix Computations},
+  Author                   = {Golub, Gene H. and Van Loan, Charles F.},
+  Publisher                = {The Johns Hopkins University Press},
+  Year                     = {1996},
+  Edition                  = {Third}
+}
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/images/9_nat.png b/01_Direct-Methods-for-Solving-Linear-Systems/Report/images/9_nat.png
new file mode 100644
index 0000000..8e93cc1
Binary files /dev/null and b/01_Direct-Methods-for-Solving-Linear-Systems/Report/images/9_nat.png differ
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/images/9_our.png b/01_Direct-Methods-for-Solving-Linear-Systems/Report/images/9_our.png
new file mode 100644
index 0000000..76b28b6
Binary files /dev/null and b/01_Direct-Methods-for-Solving-Linear-Systems/Report/images/9_our.png differ
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/images/electro.png b/01_Direct-Methods-for-Solving-Linear-Systems/Report/images/electro.png
new file mode 100644
index 0000000..a02fc5f
Binary files /dev/null and b/01_Direct-Methods-for-Solving-Linear-Systems/Report/images/electro.png differ
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/images/logo-pwr-2016.pdf b/01_Direct-Methods-for-Solving-Linear-Systems/Report/images/logo-pwr-2016.pdf
new file mode 100644
index 0000000..c919173
Binary files /dev/null and b/01_Direct-Methods-for-Solving-Linear-Systems/Report/images/logo-pwr-2016.pdf differ
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/main.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/main.tex
new file mode 100644
index 0000000..3d0e81b
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/main.tex
@@ -0,0 +1,108 @@
+\documentclass[a4paper]{article}
+\usepackage{WUST-AAE-NMaO-Report}
+\usepackage{natbib}
+
+\usepackage{amsmath}
+% \usepackage{amssymb}
+% \usepackage[polish,main=english]{babel}
+% \usepackage{titlesec}
+% \usepackage{vmargin}
+% \usepackage{graphicx}
+% \graphicspath{{images/}}
+% \usepackage{fancyhdr}
+% \usepackage{systeme}
+% \usepackage{xcolor}
+% \usepackage{indentfirst}
+% \usepackage{xspace}
+% \newcommand{\MATLAB}{\textsc{Matlab}\xspace}
+% \usepackage{listings}
+% \usepackage{matlab-prettifier}
+% \usepackage{hyperref}
+% \usepackage{multicol}
+
+
+% \newcommand{\sectionbreak}{\clearpage}
+% \setmarginsrb{3 cm}{2.5 cm}{3 cm}{2.5 cm}{1 cm}{1.5 cm}{1 cm}{1.5 cm}
+
+% Augmented matrix
+% \newenvironment{amatrix}[2]{%
+  % \left[\begin{array}{@{}*{#1}{c} | c *{#2}{c} @{}}
+% }{%
+  % \end{array}\right]
+% }
+
+% Bold letters in matrices
+% \newcommand{\matr}[1]{\mathbf{#1}}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\title{Numerical Algorithms – Report 1: Direct Methods for Solving Linear Systems}
+\author{Sergiusz Warga 230757}
+\date{\today}
+\reporttutor{dr hab. inż. Rafał Zdunek}
+
+\begin{document}
+
+\maketitle
+\tableofcontents
+\pagebreak
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Problems calculations}
+\input{problems/Problem_1}
+% \input{problems/Problem2}\\
+% \input{problems/Problem3}\\
+% \input{problems/Problem4}\\
+% \input{problems/Problem5}\\
+% \input{problems/Problem6}\\
+% \input{problems/Problem7}\\
+% \input{problems/Problem8}\\
+% \input{problems/Problem9}\\
+% \input{problems/Problem10}\\
+% \input{problems/Problem11}\\
+% \input{problems/Problem12}\\
+% \input{problems/Problem13}\\
+% \input{problems/Problem14}\\
+% \input{problems/Problem15}
+
+% \section{Algorithms}
+% \subsection{Algorithm 1 -- Gaussian elimination}\label{algorithm:1}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg1.m}
+% \subsection{Algorithm 2 -- Gaussian  elimination  with  Complete  Pivoting}\label{algorithm:2}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg2.m}
+% \subsection{Algorithm 3 -- Forward Substitution}\label{algorithm:3}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg3.m}
+% \subsection{Algorithm 4 -- Back Substitution}\label{algorithm:4}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg4.m}
+% \subsection{Algorithm 5 -- The Gauss-Jordan elimination algorithm}\label{algorithm:5}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg5.m}
+% \subsection{Algorithm 6 -- The RREF algorithm}\label{algorithm:6}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg6_RREF.m}
+% \subsection{Algorithm 7 -- The LU  factorization without pivoting}\label{algorithm:7}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg7.m}
+% \subsection{Algorithm 8 -- The LU  factorization with partial pivoting}\label{algorithm:8}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg8.m}
+% % \subsection{Algorithm 9 -- A family of the Cholesky  factorization  algorithms}\label{algorithm:9}
+% % \lstinputlisting[style=Matlab-editor]{algorithms/Alg9.m}
+% \stepcounter{subsection}
+% \subsection{Algorithm 10 -- The Cholesky factorization}\label{algorithm:10}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg10.m}
+% \subsection{Algorithm 11 -- The  QR  algorithm  by  the  Householder  transformation}\label{algorithm:11}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg11.m}
+% \subsection{Algorithm 12 -- The QR algorithm by the Givens rotations}\label{algorithm:12}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg12.m}
+% \subsection{Algorithm 13 -- The QR algorithm by the Gram-Schmidt orthogonalization}\label{algorithm:13}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg13.m}
+
+%%%%%%%%%%%%%%%%%%%
+%% BIBLIOGRAPHY %%%
+%%%%%%%%%%%%%%%%%%%
+
+\clearpage
+
+\nocite{Zdunek, GoluVanl96}
+\bibliographystyle{alpha}
+\bibliography{bibliography}
+
+\end{document}
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem1.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem1.m
new file mode 100644
index 0000000..3db2b87
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem1.m
@@ -0,0 +1,13 @@
+A = [2, -1, 0, 0;
+    -1, 2, -1, 0;
+    0, -1, 2, -1;
+    0, 0, -1, 2];
+b = [0;0;0;5];
+B = Alg5([A b])
+
+B =
+
+    1.0000         0         0         0    1.0000
+         0    1.0000         0         0    2.0000
+         0         0    1.0000         0    3.0000
+         0         0         0    1.0000    4.0000
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem10.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem10.tex
new file mode 100644
index 0000000..b288704
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem10.tex
@@ -0,0 +1,58 @@
+\subsection{Problem 10}
+Transform  the  following  matrix  to  the  RREF,  determine $rank(\mathbf{A})$ and  identify  the columns corresponding to the basic and free variables. 
+\begin{equation*}
+    \matr{A} = 
+    \begin{bmatrix}
+    1 & 2 & 2 & 3 & 1 \\
+    2 & 4 & 4 & 6 & 2 \\
+    3 & 6 & 6 & 9 & 6 \\
+    1 & 2 & 4 & 5 & 3 
+    \end{bmatrix}
+\end{equation*}
+
+Check whether the symmetric matrix $\mathbf{A}$ is positive-definite. If so, apply the Cholesky factorization. Then, compute its inverse.
+
+\subsubsection*{Solution}
+
+\begin{equation*}
+    \begin{split}
+    \begin{bmatrix}
+    1 & 2 & 2 & 3 & 1 \\
+    2 & 4 & 4 & 6 & 2 \\
+    3 & 6 & 6 & 9 & 6 \\
+    1 & 2 & 4 & 5 & 3 
+    \end{bmatrix}
+    \xrightarrow{\substack{R_2 - 2R_1 \\ R_3 - 3R_1 \\ R_4 - R_1}}
+    \begin{bmatrix}
+    1 & 2 & 2 & 3 & 1 \\
+    0 & 0 & 0 & 0 & 0 \\
+    0 & 0 & 0 & 0 & 3 \\
+    0 & 0 & 2 & 2 & 2
+    \end{bmatrix}
+    \xrightarrow{pivot}
+    \begin{bmatrix}
+    1 & 2 & 2 & 3 & 1 \\
+    0 & 0 & 2 & 2 & 2 \\
+    0 & 0 & 0 & 0 & 3 \\
+    0 & 0 & 0 & 0 & 0
+    \end{bmatrix}
+    \\
+    \xrightarrow{\substack{R_2 / 2 \\ R_3 / 3}}
+    \begin{bmatrix}
+    1 & 2 & 2 & 3 & 1 \\
+    0 & 0 & 1 & 1 & 1 \\
+    0 & 0 & 0 & 0 & 1 \\
+    0 & 0 & 0 & 0 & 0
+    \end{bmatrix}
+    \xrightarrow{\substack{R_1 - R_3 \\ R_2 - R_3}}
+    \begin{bmatrix}
+    1 & 2 & 2 & 3 & 0 \\
+    0 & 0 & 1 & 1 & 0 \\
+    0 & 0 & 0 & 0 & 1 \\
+    0 & 0 & 0 & 0 & 0
+    \end{bmatrix}
+    \end{split}
+\end{equation*}
+
+Rank of the given matrix is equal to 3.
+
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem11.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem11.tex
new file mode 100644
index 0000000..21a7f33
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem11.tex
@@ -0,0 +1,17 @@
+\subsection{Problem 11}
+
+Compute the LU factorization for
+
+\begin{equation*}
+    \matr{A} = 
+    \begin{bmatrix}
+    \phantom{-}1 & \phantom{-}3 & 3 & 2 \\
+    \phantom{-}2 & \phantom{-}6 & 9 & 5 \\
+    -1 & -3 & 3 & 0
+    \end{bmatrix}
+\end{equation*}
+
+
+
+
+Determine a set of basic variables and a set of free variables, and find a homogeneous solution to $\mathbf{Ax = 0}$. What is the rank of $\mathbf{A}$? 
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem12.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem12.tex
new file mode 100644
index 0000000..bc33872
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem12.tex
@@ -0,0 +1,14 @@
+\subsection{Problem 12}
+
+Compute the QR factorization of the matrix:
+
+\begin{equation*}
+    \matr{A} = 
+    \begin{bmatrix}
+    0 & -1 & -3  \\
+    0 & \phantom{-}0 & -2 \\
+    0 & -2 & -1 
+    \end{bmatrix}
+\end{equation*}
+
+How many flops (multiplications/divisions, additions/subtractions) are needed to perform the QR factorization with the Householder transformations and Givens rotations?
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem13.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem13.tex
new file mode 100644
index 0000000..2be42ea
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem13.tex
@@ -0,0 +1,3 @@
+\subsection{Problem 13}
+
+Let $\mathbf{Ax = b}$ , where $A = I_N \bigotimes C^TC$, the symbol $\bigotimes$ denotes the Kronecker product, $I_N \in R^{N \times N}$ is an identity matrix, $C \in R^{M \times M}$ is a random matrix with a uniform distribution, $M = 100$, and $N=50$, and $x \square N(0, I_{MN})$. Find the direct method that solves the above system  of  linear  equations  with  the  lowest  computational  cost.  Estimate  the  cost  with  a  roughly calculated number of flops and with the elapsed time.
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem14.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem14.tex
new file mode 100644
index 0000000..678c488
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem14.tex
@@ -0,0 +1,4 @@
+\subsection{Problem 14}
+
+ Solve $\mathbf{Ax=b}$, where $A \in R ^ {10 \times 10}$ is  the  Hilbert  matrix, and,  $x \square N(0, I_{10})$ with the Gaussian elimination in Matlab. Then make a small change in an entry of $\mathbf{A}$ or $\mathbf{b}$, and compare the solutions.
+  
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem15.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem15.tex
new file mode 100644
index 0000000..a1258fe
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem15.tex
@@ -0,0 +1,59 @@
+\subsection{Problem 15}
+
+Find the currents in the circuit, assuming all resistances are 10 Ohm.
+\begin{figure}[h]
+    \centering
+    \includegraphics[width=0.5\textwidth]{images/electro.png}
+    \label{fig:electro}
+\end{figure}
+
+\subsubsection*{Solution}
+
+\begin{equation*}
+    \left\{\begin{matrix}
+I_{L1}(R_1 + R_2) - I_{L3}(R_2) = -E_2\\ 
+I_{L2}(R_3 + R_4) - I_{L3}(R_3) = E_2\\ 
+-I_{L1}(R_2) - I_{L2}(R_3) + I_{L3}(R_2 + R_3) = -E_1
+\end{matrix}\right.
+\end{equation*}
+
+After substitution we have:
+
+\begin{equation*}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+    \end{amatrix} = 
+    \begin{amatrix}{3}{1}
+        20 & 0 & -10 & -10 \\
+        0 & 20 & -10 & 10 \\
+        -10 & -10 & 20 & -20
+    \end{amatrix}
+    \xrightarrow[magic]{}
+    \begin{amatrix}{3}{1}
+        1 & 0 & 0 & -1.5 \\
+        0 & 1 & 0 & -0.5 \\
+        0 & 0 & 1 & -2
+    \end{amatrix}
+\end{equation*}
+
+\begin{equation*}
+\left\{\begin{matrix}
+I_1 = -I_{L3} \\
+I_2 = -I_{L1} + I_{L2} \\
+I_3 = I_{L1} - I_{L3} \\
+I_4 = -I_{L1} \\
+I_5 = I_{L2} - I_{L3} \\
+I_6 = -I_{L2}
+\end{matrix}\right.
+\end{equation*}
+
+\begin{equation*}
+\left\{\begin{matrix}
+I_1 = 2 \\
+I_2 = 1 \\
+I_3 = 0.5 \\
+I_4 = 1.5 \\
+I_5 = 1.5 \\
+I_6 = 0.5
+\end{matrix}\right.
+\end{equation*}
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem2.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem2.m
new file mode 100644
index 0000000..e4fc915
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem2.m
@@ -0,0 +1,15 @@
+A = [1, 1, 1;
+     1, 1, 2;
+     1, 2, 2];
+b = [1;2;1];
+[P, Q, L, U] = Alg2(A);
+% PAQ = LU
+% Ly = b and Ux = y
+y = Alg3(L, P*b);     % Forward substitution
+x = Q*Alg4(U, y)      % Backward substitution
+
+x =
+
+     1
+    -1
+     1
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem2.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem2.tex
new file mode 100644
index 0000000..510b40b
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem2.tex
@@ -0,0 +1,52 @@
+\subsection{Problem 2}\label{problem:2}
+Solve the following system of linear equations using the Gaussian elimination with pivoting:
+\begin{equation*}
+    \systeme{x_1+x_2+x_3=1,x_1+x_2+2x_3=2,x_1+2x_2+2x_3=1}
+\end{equation*}
+Explain why the Gaussian elimination without pivoting does not work.
+\subsubsection*{Mathematics}
+As was stated at the end of the theoretical introduction to the previous problem, Gaussian Elimination algorithms require the diagonal entry to be different than zero (to not normalize the row dividing by zero). While in manual computations such a division would be noticed, in \MATLAB their result is \texttt{Inf} and algorithm proceeds with replacing subsequent coefficients with this value, resulting in gibberish output. With pivoting this behaviour can be avoided. The basic idea is to find a maximum absolute value in the zeroed column below the pivot and interchange their rows. This is called a partial pivoting. Searching in the whole sub-matrix spanned between pivot and $A_{n,m}$, therefore interchanging not only rows but also columns, is called a complete pivoting.
+\subsubsection*{Solution}
+\begin{equation*}
+\begin{split}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+    \end{amatrix} = 
+    \begin{amatrix}{3}{1}
+        1 & 1 & 1 & 1\\
+        1 & 1 & 2 & 2\\
+        1 & 2 & 2 & 1
+    \end{amatrix} \xrightarrow[\substack{R_2-R_1\\R_3-R_1}]{}
+    \begin{amatrix}{3}{1}
+        1 & 1 & 1 & 1\\
+        0 & 0 & 1 & 1\\
+        0 & 1 & 1 & 0
+    \end{amatrix} \xrightarrow[R_2\leftrightarrow R_3]{}
+    \begin{amatrix}{3}{1}
+        1 & 1 & 1 & 1\\
+        0 & 1 & 1 & 0\\
+        0 & 0 & 1 & 1
+    \end{amatrix}\rightarrow\\
+    \rightarrow
+    \systeme{x_1+x_2+x_3=1,x_2+x_3=0,x_3=1}
+    \rightarrow
+    \systeme{x_1+x_2+x_3=\phantom{-}1,x_2=-1,x_3=\phantom{-}1}
+    \rightarrow
+    \begin{cases}
+    x_1=\phantom{-}1\\
+    x_2=-1\\
+    x_3=\phantom{-}1
+    \end{cases}
+\end{split}
+\end{equation*}
+
+To verify manual solution the Algorithm~\ref{algorithm:2} was used. As is explained in detail in~\cite[section 3.4.2]{GoluVanl96} it transforms matrix $\matr{P}\matr{A}\matr{U}$ into a form from which lower unit triangular and upper triangular matrices may be extracted, so that:
+\begin{align*}
+   \matr{P}\matr{A}\matr{U} &= \matr{L}\matr{U} & \matr{A}\matr{x}&=\matr{b}
+\end{align*}
+therefore
+\begin{align*}
+    \matr{A}&=\matr{P^T}\matr{L}\matr{U}\matr{Q^T} & \matr{P^T}\matr{L}\matr{U}\matr{Q^T}\matr{x} &= \matr{b}\\
+    & & \matr{L}\matr{U}\matr{Q^T}\matr{x} &= \matr{P}\matr{b}
+\end{align*}
+\lstinputlisting[style=Matlab-editor]{problems/Problem2.m}
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem3.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem3.tex
new file mode 100644
index 0000000..a2cc5c7
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem3.tex
@@ -0,0 +1,140 @@
+\subsection{Problem 3}
+Solve the system:
+\begin{equation*}
+    \systeme{0.0001x_1+x_2=1,x_1+x_2=2}
+\end{equation*}
+with the Gaussian elimination with and without pivoting at the round-off error limited to 3 significant digits. Compute the condition number of the system matrix.
+\subsubsection*{Mathematics}
+Not only pivoting helps with $a_{ii}=0$ problem but, by choosing the biggest entry, ensures the stability of the LU factorization. As explained in detail in~~\cite[section 3.3]{GoluVanl96} for the factorization errors to be small, the diagonal entries of $\matr{U}$ matrix has to be small. Please note that in both solutions and condition number computations precision was limited to 3 significant digits.
+\subsubsection*{Solution without pivoting}
+\begin{equation*}
+\begin{split}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+    \end{amatrix} = 
+    \begin{amatrix}{2}{1}
+        10^{-4} & 1 & 1\\
+        1 & 1 & 2
+    \end{amatrix}
+    \xrightarrow[\substack{R_2-10^{4}R_1}]{}
+    \begin{amatrix}{2}{1}
+        10^{-4} & 1 & 1\\
+        0 & -10^{4} & -10^{4}
+    \end{amatrix}
+\end{split}
+\end{equation*}
+gives us
+\begin{equation*}
+    \begin{cases}
+    x_1=0\\
+    x_2=1\\
+    \end{cases}
+\end{equation*}
+which is obviously contradictory to $R_2$.
+
+\subsubsection*{Solution with pivoting}
+\begin{equation*}
+\begin{split}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+    \end{amatrix} = 
+    \begin{amatrix}{2}{1}
+        10^{-4} & 1 & 1\\
+        1 & 1 & 2
+    \end{amatrix}
+    \xrightarrow[R_1\leftrightarrow R_2]{}
+    \begin{amatrix}{2}{1}
+        1 & 1 & 2 \\
+        10^{-4} & 1 & 1
+    \end{amatrix}
+    \xrightarrow[\substack{R_2-10^{-4}R_1}]{}
+    \begin{amatrix}{2}{1}
+        1 & 1 & 2 \\
+        0 & 1 & 1
+    \end{amatrix}
+\end{split}
+\end{equation*}
+gives us
+\begin{equation*}
+    \begin{cases}
+    x_1=1\\
+    x_2=1
+    \end{cases}
+\end{equation*}
+\subsubsection*{Condition number}
+Condition number should describe the stability of a computed matrix. The bigger the condition number $\kappa$ the less stable the (coefficient) matrix is. It is computed as: 
+\begin{equation}
+    \kappa(\mathbf{A})=||\mathbf{A}||_2\cdot||\mathbf{A}^{-1}||_2
+\end{equation}
+where 
+\begin{equation}\label{eqn:2-norm}
+\begin{split}
+    ||\mathbf{A}||_2 = \sqrt{\lambda_{max}A^TA}
+\end{split}
+\end{equation}
+$\matr{A}^{-1}$ may be computed via Gauss-Jordan algorithm with pivoting:
+\begin{equation*}
+\begin{split}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{I}
+    \end{amatrix} = 
+    \begin{amatrix}{2}{2}
+        \frac{1}{1000} & 1 & 1 & 0\\
+        1 & 1 & 0 & 1
+    \end{amatrix}
+    \xrightarrow[R_1\leftrightarrow R_2]{}
+    \begin{amatrix}{2}{2}
+        1 & 1 & 0 & 1\\
+        \frac{1}{1000} & 1 & 1 & 0
+    \end{amatrix}
+    \xrightarrow[\substack{R_2-\frac{1}{1000}R_2}]{}\\
+    \xrightarrow[\substack{R_2-\frac{1}{1000}R_2}]{}
+    \begin{amatrix}{2}{2}
+        1 & 1 & 0 & 1\\
+        0 & 1 & 1 & 0
+    \end{amatrix}
+    \xrightarrow[\substack{R_1-R_2}]{}
+    \begin{amatrix}{2}{2}
+        1 & 0 & -1 & 1\\
+        0 & 1 & 1 & 0
+    \end{amatrix} = 
+    \begin{amatrix}{1}{1}
+        \matr{I} & \matr{A^{-1}}
+    \end{amatrix}
+\end{split}
+\end{equation*}
+Now, using the equation~\ref{eqn:2-norm}
+
+\parbox{0.5\textwidth}{
+\begin{gather*}
+    \matr{A}^T\matr{A} = 
+    \begin{bmatrix}
+        1 & 1 \\
+        1 & 2
+    \end{bmatrix} \\
+    \det\begin{pmatrix}
+        1-\lambda & 1 \\
+        1 & 2-\lambda
+    \end{pmatrix} = 0 \\
+    \lambda=\frac{3\pm\sqrt{5}}{2}
+    \end{gather*}
+}
+\parbox{0.5\textwidth}{
+\begin{gather*}
+     \matr{A}^{-1^T}\matr{A}^{-1} = 
+    \begin{bmatrix}
+        2 & -1 \\
+        -1 & 1
+    \end{bmatrix} \\
+    \det\begin{pmatrix}
+        2-\lambda & -1 \\
+        -1 & 1-\lambda
+    \end{pmatrix} = 0 \\
+    \lambda=\frac{3\pm\sqrt{5}}{2}
+\end{gather*}
+}
+we finally obtain
+\begin{equation*}
+    \kappa(\mathbf{A})=\frac{3\pm\sqrt{5}}{2}\approx2.618
+\end{equation*}
+which is a very low number indicating, that the matrix $\matr{A}$ is computationally stable.
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem4.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem4.m
new file mode 100644
index 0000000..6062280
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem4.m
@@ -0,0 +1,23 @@
+A = [0.835, 0.667;
+     0.333, 0.266];
+b = [0.168; 0.067];
+bp = [0.168; 0.066];
+B = Alg5([A b])
+Bp = Alg5([A bp])
+kappa = cond(A)
+
+B =
+
+    1.0000         0    1.0000
+         0    1.0000   -1.0000
+
+
+Bp =
+
+    1.0000         0 -666.0000
+         0    1.0000  834.0000
+
+
+kappa =
+
+   1.3238e+06
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem4.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem4.tex
new file mode 100644
index 0000000..eb239b2
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem4.tex
@@ -0,0 +1,37 @@
+\subsection{Problem 4}
+Solve the system of linear equations:
+\begin{equation*}
+    \systeme{0.835x_1 + 0.667x_2 = 0.168,0.333x_1 + 0.266x_2 = 0.067}
+\end{equation*}
+
+Then slightly perturb $b_2$ from $0.067$ to $0.066$ and compute the solution to the perturbed system. Explain the change in the solution by computing the condition number of the system matrix.
+\subsubsection*{Solution}
+\begin{equation*}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+    \end{amatrix} = 
+    \begin{amatrix}{2}{1}
+        0.835 & 0.667 & 0.168\\
+        0.333 & 0.266 & 0.067
+    \end{amatrix} \xrightarrow[R_2-0.333R_1]{\frac{R_1}{0.835}}
+    \begin{amatrix}{2}{1}
+        1 & 0.7988 & 0.2012\\
+        0 & 0 & 0
+    \end{amatrix}
+\end{equation*}
+
+which yields infinite many solutions. Perturbing $b_2$ we obtain:
+\begin{equation*}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b'}
+    \end{amatrix} \xrightarrow[R_2-0.333R_1]{\frac{R_1}{0.835}}
+    \begin{amatrix}{2}{1}
+        1 & 0.7988 & 0.2012\\
+        0 & 0 & -0.001
+    \end{amatrix}
+\end{equation*}
+
+which yields none solutions. Performing the same calculations with Algorithm~\ref{algorithm:5}:
+\lstinputlisting[style=Matlab-editor]{problems/Problem4.m}
+
+With enough precision \MATLAB presented two valid solutions. Both of them differed by a factor of hundreds, which may be explained by a huge coefficient number (way over a million).
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem5.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem5.m
new file mode 100644
index 0000000..d2bddf6
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem5.m
@@ -0,0 +1,17 @@
+A = [2, 1, 2;
+     1, 2, 3;
+     4, 1, 2];
+[P, Q, L, U] = Alg2(A);
+I = P*eye(3);
+% Ly = b and Ux = y
+y = Alg3(L, I);     % Forward substitution
+x = Q*Alg4(U, y);   % Backward substitution
+inv(A) - x
+
+ans =
+
+   1.0e-15 *
+
+   -0.2220    0.0555    0.1110
+    0.8882   -0.4441   -0.6661
+         0    0.1110         0
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem5.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem5.tex
new file mode 100644
index 0000000..b28c2c2
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem5.tex
@@ -0,0 +1,17 @@
+\subsection{Problem 5}
+Find $\mathbf{A^{-1}}$ to:
+
+\begin{equation*}
+    \matr{A} = 
+    \begin{bmatrix}
+        2 & 1 & 2 \\
+        1 & 2 & 3 \\
+        4 & 1 & 2 
+    \end{bmatrix}
+\end{equation*}
+
+solving the system $\mathbf{AX=I_{3}}$.
+
+\subsubsection*{Solution}
+\lstinputlisting[style=Matlab-editor]{problems/Problem5.m}
+which judges well our algorithms.
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem6.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem6.m
new file mode 100644
index 0000000..7cb996d
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem6.m
@@ -0,0 +1,26 @@
+A = [1 2 3 4;
+    -1 1 2 1;
+     0 2 1 3;
+     0 0 1 1];
+b = [1; 1; 1; 1];
+% PA = LU
+[L, U, P] = Alg8(A);
+det(A) - prod(diag(P.'*U))
+% Ly = b and Ux = y
+y = Alg3(L, P*b);   % Forward substitution
+x = Alg4(U, y);   % Backward substitution
+x - A\b
+
+ans =
+
+     0
+
+
+ans =
+
+   1.0e-15 *
+
+    0.4441
+         0
+         0
+         0
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem6.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem6.tex
new file mode 100644
index 0000000..b8ea8bc
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem6.tex
@@ -0,0 +1,26 @@
+\subsection{Problem 6}
+
+Apply the LU factorization to the matrix
+\begin{equation*}
+    \matr{A} = 
+    \begin{bmatrix}
+        \phantom{-}1 & 2 & 3 & 4 \\
+        -1 & 1 & 2 & 1 \\
+        \phantom{-}0 & 2 & 1 & 3 \\
+        \phantom{-}0 & 0 & 1 & 1
+    \end{bmatrix}
+\end{equation*}
+Then calculate $\det(\matr{A})$ using the matrix $\matr{U}$. Finally solve $\matr{A}\matr{x}=\matr{b}$ for $\matr{b}=[1\dots1]^T$.
+\subsubsection*{Mathematics}
+The LU Factorization decomposes matrix $\matr{A}$ into an upper triangular matrix $\matr{U}$ and unit lower triangular matrix $\matr{L}$, so that
+\begin{equation*}
+    \matr{A} = \matr{L}\matr{U}
+\end{equation*}
+Using the property of a determinant and a fact, that $\det(\matr{L}) = 1$ one can calculate $\det(\matr{A})$ as
+\begin{equation*}
+    \det(\matr{A}) = \det(\matr{L}\matr{U}) = \det(\matr{L}) \cdot \det(\matr{U}) = \det(\matr{U}) = u_{11}\cdots u_{nn}
+\end{equation*}
+\subsubsection*{Solution}
+A glance at the matrix $\matr{A}$ tells us that it is not strictly diagonally dominant, so a pivoting algorithm should be used.
+%\footnote{$\matr{A}\in\mathbb{R}^{n\times n}$ is \textit{strictly diagonally dominant} if $|a_{ii}|>\sum_{j=1 j i}^n|a_{ij}|$}
+\lstinputlisting[style=Matlab-editor]{problems/Problem6.m}
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem7.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem7.m
new file mode 100644
index 0000000..83bae51
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem7.m
@@ -0,0 +1,24 @@
+A = hilb(5);
+[L, U] = Alg7(A);
+cond(A)
+L*U - A
+prod(diag(U)) - det(A)
+
+ans =
+
+   4.7661e+05
+   
+ans =
+
+   1.0e-16 *
+
+         0         0         0         0         0
+         0         0         0         0         0
+         0         0         0         0         0
+         0         0         0         0   -0.1388
+         0         0         0   -0.1388         0
+
+
+ans =
+
+   1.6620e-23
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem7.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem7.tex
new file mode 100644
index 0000000..08c50e5
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem7.tex
@@ -0,0 +1,21 @@
+\subsection{Problem 7}
+
+Let $\mathbf{A} = \left[a_{ij}\right] \in\mathbb{R}^{N\times N}$ with, $a_{ij} = \frac{1}{i + j - 1}$ (Hilbert  matrix). For  $N=5$, perform  the LU factorization of the matrix $\mathbf{A}$. Then, compute det($\mathbf{A}$).
+\subsubsection*{Mathematics}
+Hilbert matrices are an example of ill-conditioned matrices, which -- with their high $\kappa$ -- makes the numerical computations highly unstable. For our  $\mathbf{H} = \left[h_{ij}\right] \in\mathbb{R}^{5\times 5}$:
+
+\begin{equation*}
+    \matr{H} = 
+    \begin{bmatrix}
+               1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \\
+     \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} \\
+     \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} \\
+     \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} \\
+     \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} 
+    \end{bmatrix}
+\end{equation*}
+
+$\kappa(\matr{H})\approx4.7\cdot10^5$ which is a huge value similar to this in the Problem 4.
+\subsubsection*{Solution}
+\lstinputlisting[style=Matlab-editor]{problems/Problem7.m}
+The above results may look promising, but a determinant calculated both by \MATLAB and our algorithms is of the order of $4\cdot10^{-12}$, which for numerical computations is highly not useful.
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem8.m b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem8.m
new file mode 100644
index 0000000..1128d60
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem8.m
@@ -0,0 +1,32 @@
+A = [ 1 -1  0  0;
+     -1  2 -1  0;
+      0 -1  2 -1;
+      0  0 -1  2];
+if all(eig(A) > 0)
+    C = Alg10(A)
+    C*C.' - A
+    inv(C)
+end
+
+C =
+
+     1     0     0     0
+    -1     1     0     0
+     0    -1     1     0
+     0     0    -1     1
+
+
+ans =
+
+     0     0     0     0
+     0     0     0     0
+     0     0     0     0
+     0     0     0     0
+
+
+ans =
+
+     1     0     0     0
+     1     1     0     0
+     1     1     1     0
+     1     1     1     1
\ No newline at end of file
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem8.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem8.tex
new file mode 100644
index 0000000..194a377
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem8.tex
@@ -0,0 +1,16 @@
+\subsection{Problem 8}
+Let
+\begin{equation*}
+    \matr{A} = 
+    \begin{bmatrix}
+        \phantom{-}1 & -1 & \phantom{-}0 & \phantom{-}0 \\
+        -1 & \phantom{-}2 & -1 & \phantom{-}0 \\
+        \phantom{-}0 & -1 & \phantom{-}2 & -1 \\
+        \phantom{-}0 & \phantom{-}0 & -1 & \phantom{-}2
+    \end{bmatrix}
+\end{equation*}
+
+Check whether the symmetric matrix $\mathbf{A}$ is positive-definite. If so, apply the Cholesky factorization. Then, compute its inverse.
+
+\lstinputlisting[style=Matlab-editor]{problems/Problem8.m}
+
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem9.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem9.tex
new file mode 100644
index 0000000..03b9105
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem9.tex
@@ -0,0 +1,3 @@
+\subsection{Problem 9}
+Let  $\mathbf{A}$  be  the  Pascal  matrix  of  order  100  (\textit{pascal}  function  in Matlab).  Check  whether the matrix is positive definite. If so, apply the Cholesky factorization and give interpretation of the obtained factor.
+
diff --git a/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem_1.tex b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem_1.tex
new file mode 100644
index 0000000..a5f352c
--- /dev/null
+++ b/01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem_1.tex
@@ -0,0 +1,65 @@
+\subsection{Problem 1}
+
+Solve the system and find the pivots when:
+\begin{equation*}
+    \systeme{2u-v=0,-u+2v-w=0,-v+2w-z=0,-w+2z=5}
+\end{equation*}
+
+\subsubsection*{Mathematics}
+The solution to a linear system like the one presented above may be obtained via Gauss-Jordan Elimination. An augmented matrix
+$\matr{B} = \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+\end{amatrix}$, where $\matr{A}$ is a coefficient matrix, and $\matr{b}$ is a column vector of constant terms, is assembled. Then, the algorithm implemented in \ref{algorithm:5} is performed. Note, that none of the $a_{ii}$ can be equal to zero for this algorithm to work properly\footnote{Explanation of this property is presented in \ref{problem:2}}.
+
+\subsubsection*{Solution}
+
+Manual calculations were performed. In the following transformations a row normalization is presented above the arrow and then a zeroing is below it.
+
+\begin{equation*}
+\begin{split}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+    \end{amatrix} = 
+    \begin{amatrix}{4}{1}
+        \phantom{-}2 & -1 & \phantom{-}0 & \phantom{-}0 & 0\\
+        -1 & \phantom{-}2 & -1 & \phantom{-}0 & 0\\
+        \phantom{-}0 & -1 & \phantom{-}2 & -1 & 0\\
+        \phantom{-}0 & \phantom{-}0 & -1 & \phantom{-}2 & 5
+    \end{amatrix} \xrightarrow[R_2 + R_1]{\frac{1}{2}R_1} 
+    \begin{amatrix}{4}{1}
+        \phantom{-}1 & -\frac{1}{2} & \phantom{-}0 & \phantom{-}0 & 0\\
+        \phantom{-}0 & \phantom{-}\frac{3}{2} & -1 & \phantom{-}0 & 0\\
+        \phantom{-}0 & -1 & \phantom{-}2 & -1 & 0\\
+        \phantom{-}0 & \phantom{-}0 & -1 & \phantom{-}2 & 5
+    \end{amatrix} \xrightarrow[\substack{R_1+\frac{1}{2}R_2\\R_3+R_2}]{\frac{2}{3}R_2} \\
+    \xrightarrow[R_3+\frac{2}{3}R_2]{\frac{2}{3}R_2}
+    \begin{amatrix}{4}{1}
+        \phantom{-}1 & \phantom{-}0 & -\frac{1}{3} & \phantom{-}0 & 0\\
+        \phantom{-}0 & \phantom{-}1 & -\frac{2}{3} & \phantom{-}0 & 0\\
+        \phantom{-}0 & \phantom{-}0 & \phantom{-}\frac{4}{3} & -1 & 0\\
+        \phantom{-}0 & \phantom{-}0 & -1 & \phantom{-}2 & 5
+    \end{amatrix} \xrightarrow[\substack{R_1+\frac{1}{3}R_3\\R_2+\frac{2}{3}R_3\\R_4+R_3}]{\frac{3}{4}R_3}
+    \begin{amatrix}{4}{1}
+        \phantom{-}1 & \phantom{-}0 & \phantom{-}0 & -\frac{1}{4} & 0\\
+        \phantom{-}0 & \phantom{-}1 & \phantom{-}0 & -\frac{1}{2} & 0\\
+        \phantom{-}0 & \phantom{-}0 & \phantom{-}1 & -\frac{3}{4} & 0\\
+        \phantom{-}0 & \phantom{-}0 & \phantom{-}0 & -\frac{4}{3} & 5
+    \end{amatrix} \xrightarrow[\substack{R_1+\frac{1}{4}R_4\\R_2+\frac{1}{2}R_4\\R_3+\frac{3}{4}R_4}]{-\frac{3}{4}R_4}\\
+    \xrightarrow[\substack{R_1+\frac{1}{4}R_4\\R_2+\frac{1}{2}R_4\\R_3+\frac{3}{4}R_4}]{-\frac{3}{4}R_4}
+    \begin{amatrix}{4}{1}
+        \phantom{-}1 & \phantom{-}0 & \phantom{-}0 & \phantom{-}0 & 1\\
+        \phantom{-}0 & \phantom{-}1 & \phantom{-}0 & \phantom{-}0 & 2\\
+        \phantom{-}0 & \phantom{-}0 & \phantom{-}1 & \phantom{-}0 & 3\\
+        \phantom{-}0 & \phantom{-}0 & \phantom{-}0 & \phantom{-}1 & 4
+    \end{amatrix} \rightarrow
+    \begin{cases}
+    u=1\\
+    v=2\\
+    w=3\\
+    z=4
+    \end{cases}
+\end{split}
+\end{equation*}
+
+Above result was verified with the Algorithm~\ref{algorithm:5}:
+\lstinputlisting[style=Matlab-editor]{problems/Problem1.m}
diff --git a/02_Iterative-Linear-Solvers/Report/main.tex b/02_Iterative-Linear-Solvers/Report/main.tex
new file mode 100644
index 0000000..ee6dbab
--- /dev/null
+++ b/02_Iterative-Linear-Solvers/Report/main.tex
@@ -0,0 +1,15 @@
+\documentclass[a4paper]{article}
+\usepackage{WUST-AAE-NMaO-Report}
+
+\title{Numerical Algorithms – Report 1: Iterative Linear Solvers}
+\author{Sergiusz Warga 230757}
+\date{\today}
+\reporttutor{dr hab. inż. Rafał Zdunek}
+
+\begin{document}
+
+\maketitle
+\tableofcontents
+\pagebreak
+
+\end{document}
diff --git a/Direct Methods for Solving Linear Systems/README.md b/Direct Methods for Solving Linear Systems/README.md
deleted file mode 100644
index 6119207..0000000
--- a/Direct Methods for Solving Linear Systems/README.md	
+++ /dev/null
@@ -1,8 +0,0 @@
-# Direct Methods for Solving Linear Systems
-
-## Bibliography
-
-References in the code comments are represented by `[number]`:
-
-  - 1
-  - 2 – Golub, van Loan, Matrix Computations, 3rd edition
\ No newline at end of file
diff --git a/Dockerfile b/Dockerfile
new file mode 100644
index 0000000..3fcd7d1
--- /dev/null
+++ b/Dockerfile
@@ -0,0 +1,6 @@
+FROM alpine:3.16.0
+
+RUN apk update && \
+    apk add --no-cache \
+    texlive-full=20220403.62885-r2 \
+    make=4.3-r0
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..5c93f45
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,13 @@
+            DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
+                    Version 2, December 2004
+
+ Copyright (C) 2004 Sam Hocevar <sam@hocevar.net>
+
+ Everyone is permitted to copy and distribute verbatim or modified
+ copies of this license document, and changing it is allowed as long
+ as the name is changed.
+
+            DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
+   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+  0. You just DO WHAT THE FUCK YOU WANT TO.
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..349e687
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,63 @@
+BUILDDIR = `pwd`/build
+DOCKER_IMAGE = latex-builder
+LATEX_OPTIONS = -pdf
+#  -interaction=nonstopmode
+PACKAGE_NAME = WUST-AAE-NMaO-Report
+
+SOURCE_DIRS = $(shell ls -d */ | grep -v "^build")
+SOURCE_TEXS = $(SOURCE_DIRS:=Report/main.tex)
+PDFBUILDS = $(addprefix $(BUILDDIR)/, $(SOURCE_DIRS:/=.pdf))
+
+
+.DEFAULT_GOAL := docker
+
+
+.PHONY : clean
+clean :
+	-rm -r $(BUILDDIR) $(PACKAGE_NAME).sty
+
+
+.PHONY : all
+all : $(PDFBUILDS)
+
+
+.PHONY : sty
+sty : $(PACKAGE_NAME).sty
+$(PACKAGE_NAME).sty : $(PACKAGE_NAME).ins $(PACKAGE_NAME).dtx
+	mkdir -p $(BUILDDIR)
+	yes | latex -output-directory=$(BUILDDIR) $(PACKAGE_NAME).ins
+	cp $(BUILDDIR)/$@ $@
+
+
+$(PDFBUILDS) : $(BUILDDIR)/%.pdf: $(PACKAGE_NAME).sty %/Report/main.tex
+	mkdir -p $(BUILDDIR)
+	cp logo-pwr-2016.pdf $(BUILDDIR)/
+	latexmk $(LATEX_OPTIONS) \
+		-cd \
+		-jobname=$* \
+		-output-directory=$(BUILDDIR) \
+		$*/Report/main
+
+
+.PHONY : hadolint
+hadolint : Dockerfile
+	docker run --rm --interactive hadolint/hadolint < Dockerfile
+
+
+.PHONY : image
+image : Dockerfile
+	docker build \
+		--tag $(DOCKER_IMAGE) \
+		.
+
+
+.PHONY : docker
+docker : image ## Compile the project via the latex-builder docker image
+	docker run \
+		--rm \
+		--interactive \
+		--workdir /data \
+		--volume `pwd`:/data \
+		--name=$(DOCKER_IMAGE) \
+		$(DOCKER_IMAGE) \
+		sh -c "make all"
diff --git a/README.md b/README.md
index e69de29..fa24a10 100644
--- a/README.md
+++ b/README.md
@@ -0,0 +1 @@
+# Numerical Methods and Optimization
diff --git a/WUST-AAE-NMaO-Report.dtx b/WUST-AAE-NMaO-Report.dtx
new file mode 100644
index 0000000..a6dd6af
--- /dev/null
+++ b/WUST-AAE-NMaO-Report.dtx
@@ -0,0 +1,145 @@
+% \iffalse meta-comment
+%
+% DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
+%
+% Author: Sergiusz Warga <sergiusz.warga@gmail.com>
+%
+% \fi
+%
+% \iffalse
+%<*driver>
+\ProvidesFile{WUST-AAE-NMaO-Report.dtx}
+%</driver>
+%<package>\NeedsTeXFormat{LaTeX2e}[2005/12/01]
+%<package>\ProvidesPackage{WUST-AAE-NMaO-Report}
+%<*package>
+  [2023/01/06 v0.1.0 WUST AAE NMaO Report package]
+%</package>
+%
+%<*driver>
+\documentclass{ltxdoc}
+\usepackage[
+  hidelinks,
+  bookmarks,
+  bookmarksopen,
+  bookmarksopenlevel=2,
+  pdftitle={WUST AAE NMaO Report LaTeX package},
+  pdfauthor={Sergiusz Warga}
+]{hyperref}
+\EnableCrossrefs
+\CodelineIndex
+\RecordChanges
+\setlength\hfuzz{15pt}      % don’t show so many
+\hbadness=7000              % over- and underfull box warnings
+\begin{document}
+  \DocInput{WUST-AAE-NMaO-Report.dtx}
+\end{document}
+%</driver>
+% \fi
+%
+% \CheckSum{0}
+%
+% \changes{v0.1.0}{2023/01/06}{Initial version}
+%
+% \GetFileInfo{WUST-AAE-NMaO-Report.dtx}
+%
+% \DoNotIndex{\@date, \@title}
+% \DoNotIndex{\\}
+% \DoNotIndex{\begin}
+% \DoNotIndex{\DeclareRobustCommand, \dimexpr, \doublerulesep}
+% \DoNotIndex{\end}
+% \DoNotIndex{\fbox, \fboxrule, \fboxsep}
+% \DoNotIndex{\gdef}
+% \DoNotIndex{\hoffset, \Huge}
+% \DoNotIndex{\Large}
+% \DoNotIndex{\newlength, \newpage, \null}
+% \DoNotIndex{\oddsidemargin}
+% \DoNotIndex{\paperwidth}
+% \DoNotIndex{\renewcommand}
+% \DoNotIndex{\setlength}
+% \DoNotIndex{\textwidth}
+% \DoNotIndex{\uppercase}
+% \DoNotIndex{\vfill}
+%
+% \title{WUST AAE NMaO Report package\thanks
+%   {This document corresponds to \texttt{WUST-AAE-NMaO-Report}~\fileversion,
+%       dated \filedate.}}\author{Sergiusz Warga\texttt{sergiusz.warga@gmail.com}}
+%
+% \maketitle
+%
+% \StopEventually{\PrintChanges \PrintIndex}
+%
+% \section{Implementation}
+%
+% \subsection{Packages}
+%    \begin{macrocode}
+\RequirePackage{graphicx}
+\RequirePackage{xspace}
+%    \end{macrocode}
+%
+% \subsection{Dimentions}
+%
+% Text area.
+%
+%    \begin{macrocode}
+\setlength{\textwidth}{140mm}
+%    \end{macrocode}
+%
+% Center the text area.
+%
+%    \begin{macrocode}
+\setlength{\hoffset}{
+  \dimexpr
+  (\paperwidth-\textwidth)/2-\oddsidemargin-1in
+}
+
+\DeclareRobustCommand*{\reportgroup}[1]{\gdef\@reportgroup{#1}}
+\DeclareRobustCommand*{\reporttutor}[1]{\gdef\@reporttutor{#1}}
+\DeclareRobustCommand*{\reportclassdate}[1]{\gdef\@reportclassdate{#1}}
+\DeclareRobustCommand*{\reportlabassistant}[1]{\gdef\@reportlabassistant{#1}}
+
+\reportgroup{}
+\reporttutor{}
+\reportclassdate{}
+\reportlabassistant{}
+
+%    \end{macrocode}
+% \DescribeMacro{\maketitle}
+%    \begin{macrocode}
+\renewcommand{\maketitle}{
+\begin{titlepage}
+  \thispagestyle{empty}
+  \begin{center}
+    \includegraphics{logo-pwr-2016.pdf}
+    \vspace{\baselineskip}
+    \rule{\linewidth}{0.2 mm}\\[0.4 cm]
+    {\huge \bfseries \@title}\\
+    \rule{\linewidth}{0.2 mm}\\[1.5 cm]
+  \vfill
+  \large
+  \textsc{\@author}\\
+  \vspace{\baselineskip}
+  \textsc{Tutor: \@reporttutor}\\
+  \textsc{Laboratory group: Monday 13:15}\\
+
+  \end{center}
+\end{titlepage}
+}
+%    \end{macrocode}
+% MATLAB command.
+%    \begin{macrocode}
+\newcommand{\MATLAB}{\textsc{Matlab}\xspace}
+%    \end{macrocode}
+% Bold letters in matrices
+%    \begin{macrocode}
+\newcommand{\matr}[1]{\mathbf{#1}}
+%    \end{macrocode}
+% Augmented matrix
+%    \begin{macrocode}
+\newenvironment{amatrix}[2]{%
+  \left[\begin{array}{@{}*{#1}{c} | c *{#2}{c} @{}}
+}{%
+  \end{array}\right]
+}
+%    \end{macrocode}
+% \Finale
diff --git a/WUST-AAE-NMaO-Report.ins b/WUST-AAE-NMaO-Report.ins
new file mode 100644
index 0000000..f94d3f6
--- /dev/null
+++ b/WUST-AAE-NMaO-Report.ins
@@ -0,0 +1,38 @@
+%%
+%% DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
+%%
+%% Author: Sergiusz Warga <sergiusz.warga@gmail.com>
+%%
+\input docstrip.tex
+\keepsilent % Turn off DocStrip activity reports
+
+\usedir{tex/latex/WUST-AAE-NMaO-Report}
+
+\preamble
+
+This is a generated file.
+
+To what the fuck you want with it.
+
+\endpreamble
+
+\generate{
+\file{WUST-AAE-NMaO-Report.sty}
+{\from{WUST-AAE-NMaO-Report.dtx}{package}}
+}
+
+\obeyspaces
+\Msg{****************************************************}
+\Msg{*                                                  *}
+\Msg{* To finish the installation you have to move the  *}
+\Msg{* following file into a directory searched by TeX: *}
+\Msg{*                                                  *}
+\Msg{*     WUST-AAE-NMaO-Report.sty                     *}
+\Msg{*                                                  *}
+\Msg{* To produce the documentation run the file        *}
+\Msg{* WUST-AAE-NMaO-Report.dtx through LaTeX.          *}
+\Msg{*                                                  *}
+\Msg{* Happy TeXing!                                    *}
+\Msg{*                                                  *}
+\Msg{****************************************************}
+\endbatchfile
diff --git a/WUST-AAE-NMaO-Report.sty b/WUST-AAE-NMaO-Report.sty
new file mode 100644
index 0000000..caa2211
--- /dev/null
+++ b/WUST-AAE-NMaO-Report.sty
@@ -0,0 +1,62 @@
+%%
+%% This is file `WUST-AAE-NMaO-Report.sty',
+%% generated with the docstrip utility.
+%%
+%% The original source files were:
+%%
+%% WUST-AAE-NMaO-Report.dtx  (with options: `package')
+%% 
+%% This is a generated file.
+%% 
+%% To what the fuck you want with it.
+%% 
+\NeedsTeXFormat{LaTeX2e}[2005/12/01]
+\ProvidesPackage{WUST-AAE-NMaO-Report}
+  [2023/01/06 v0.1.0 WUST AAE NMaO Report package]
+\RequirePackage{graphicx}
+\RequirePackage{xspace}
+\setlength{\textwidth}{140mm}
+\setlength{\hoffset}{
+  \dimexpr
+  (\paperwidth-\textwidth)/2-\oddsidemargin-1in
+}
+
+\DeclareRobustCommand*{\reportgroup}[1]{\gdef\@reportgroup{#1}}
+\DeclareRobustCommand*{\reporttutor}[1]{\gdef\@reporttutor{#1}}
+\DeclareRobustCommand*{\reportclassdate}[1]{\gdef\@reportclassdate{#1}}
+\DeclareRobustCommand*{\reportlabassistant}[1]{\gdef\@reportlabassistant{#1}}
+
+\reportgroup{}
+\reporttutor{}
+\reportclassdate{}
+\reportlabassistant{}
+
+\renewcommand{\maketitle}{
+\begin{titlepage}
+  \thispagestyle{empty}
+  \begin{center}
+    \includegraphics{logo-pwr-2016.pdf}
+    \vspace{\baselineskip}
+    \rule{\linewidth}{0.2 mm}\\[0.4 cm]
+    {\huge \bfseries \@title}\\
+    \rule{\linewidth}{0.2 mm}\\[1.5 cm]
+  \vfill
+  \large
+  \textsc{\@author}\\
+  \vspace{\baselineskip}
+  \textsc{Tutor: \@reporttutor}\\
+  \textsc{Laboratory group: Monday 13:15}\\
+
+  \end{center}
+\end{titlepage}
+}
+\newcommand{\MATLAB}{\textsc{Matlab}\xspace}
+\newcommand{\matr}[1]{\mathbf{#1}}
+\newenvironment{amatrix}[2]{%
+  \left[\begin{array}{@{}*{#1}{c} | c *{#2}{c} @{}}
+}{%
+  \end{array}\right]
+}
+\endinput
+%%
+%% End of file `WUST-AAE-NMaO-Report.sty'.
diff --git a/logo-pwr-2016.pdf b/logo-pwr-2016.pdf
new file mode 100644
index 0000000..c919173
Binary files /dev/null and b/logo-pwr-2016.pdf differ