From 71664211411474065bec4af3223e78eb751ed81f Mon Sep 17 00:00:00 2001
From: Cristiano Nunes <cfgnunes@gmail.com>
Date: Tue, 1 Mar 2022 22:46:27 -0300
Subject: [PATCH] Code refactor

---
 README.md                                     |   2 +-
 .../derivative_backward_difference.m          |  21 ++-
 differentiation/derivative_five_point.m       |  21 ++-
 differentiation/derivative_three_point.m      |  21 ++-
 integration/composite2_simpson.m              |  16 +-
 integration/composite2_trapezoidal.m          |  14 +-
 integration/composite_simpson.m               |  20 ++-
 integration/composite_trapezoidal.m           |  18 ++-
 interpolation/lagrange.m                      |  22 ++-
 interpolation/neville.m                       |  22 +--
 linear_systems/backward_substitution.m        |  16 +-
 linear_systems/forward_substitution.m         |  16 +-
 linear_systems/gauss_elimination_pp.m         |  23 ++-
 linear_systems_iterative/gauss_seidel.m       |  26 ++--
 linear_systems_iterative/jacobi.m             |  26 ++--
 main.m                                        | 140 +++++++++---------
 ode/euler.m                                   |  25 ++--
 ode/rk4.m                                     |  25 ++--
 ode/rk4_system.m                              |  25 ++--
 ode/taylor2.m                                 |  27 ++--
 ode/taylor4.m                                 |  33 +++--
 polynomials/briot_ruffini.m                   |  19 ++-
 polynomials/newton_divided_difference.m       |  27 +++-
 solutions/bisection.m                         |  25 ++--
 solutions/newton.m                            |  28 ++--
 solutions/secant.m                            |  29 ++--
 26 files changed, 407 insertions(+), 280 deletions(-)

diff --git a/README.md b/README.md
index d4d3ba8..7076678 100644
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
 
 Numerical methods implementation in MATLAB.
 
-For the "Python version", see [this repository](https://github.com/cfgnunes/numerical-methods-python).
+For the implementation in Python, see [this repository](https://github.com/cfgnunes/numerical-methods-python).
 
 ## Getting Started
 
diff --git a/differentiation/derivative_backward_difference.m b/differentiation/derivative_backward_difference.m
index e00bc20..4d47890 100644
--- a/differentiation/derivative_backward_difference.m
+++ b/differentiation/derivative_backward_difference.m
@@ -1,11 +1,14 @@
 function [dy] = derivative_backward_difference(x, y)
-% Calculate the first derivative
-% All values in 'x' must be equally spaced
-% Inputs:
-%         x: Array containing x values
-%         y: Array containing y values
-% Outputs:
-%        dy: Array containing the first derivative values
+    % Calculate the first derivative.
+    %
+    % All values in 'x' must be equally spaced.
+    %
+    % Args:
+    %     x: array containing x values.
+    %     y: array containing y values.
+    %
+    % Returns:
+    %     dy: array containing the first derivative values.
 
     x_size = size(x, 2);
     y_size = size(y, 2);
@@ -22,7 +25,9 @@
 
     n = x_size;
     dy = zeros(1, n);
+
     for i = 1:n
+
         if i == n
             h = x(i) - x(i - 1);
             dy(i) = dy_difference(- h, y(i), y(i - 1));
@@ -30,5 +35,7 @@
             h = x(i + 1) - x(i);
             dy(i) = dy_difference(h, y(i), y(i + 1));
         end
+
     end
+
 end
diff --git a/differentiation/derivative_five_point.m b/differentiation/derivative_five_point.m
index 2ed8e32..2791864 100644
--- a/differentiation/derivative_five_point.m
+++ b/differentiation/derivative_five_point.m
@@ -1,11 +1,14 @@
 function [dy] = derivative_five_point(x, y)
-% Calculate the first derivative
-% All values in 'x' must be equally spaced
-% Inputs:
-%         x: Array containing x values
-%         y: Array containing y values
-% Outputs:
-%        dy: Array containing the first derivative values
+    % Calculate the first derivative.
+    %
+    % All values in 'x' must be equally spaced.
+    %
+    % Args:
+    %     x: array containing x values.
+    %     y: array containing y values.
+    %
+    % Returns:
+    %     dy: array containing the first derivative values.
 
     x_size = size(x, 2);
     y_size = size(y, 2);
@@ -24,7 +27,9 @@
     h = x(2) - x(1);
     n = x_size;
     dy = zeros(1, n);
+
     for i = 1:n
+
         if i == 1 || i == 2
             dy(i) = dy_end(h, y(i), y(i + 1), y(i + 2), y(i + 3), y(i + 4));
         elseif i == n || i == n - 1
@@ -32,5 +37,7 @@
         else
             dy(i) = dy_mid(h, y(i - 2), y(i - 1), y(i + 1), y(i + 2));
         end
+
     end
+
 end
diff --git a/differentiation/derivative_three_point.m b/differentiation/derivative_three_point.m
index 325c5d8..976691f 100644
--- a/differentiation/derivative_three_point.m
+++ b/differentiation/derivative_three_point.m
@@ -1,11 +1,14 @@
 function [dy] = derivative_three_point(x, y)
-% Calculate the first derivative
-% All values in 'x' must be equally spaced
-% Inputs:
-%         x: Array containing x values
-%         y: Array containing y values
-% Outputs:
-%        dy: Array containing the first derivative values
+    % Calculate the first derivative.
+    %
+    % All values in 'x' must be equally spaced.
+    %
+    % Args:
+    %     x: array containing x values.
+    %     y: array containing y values.
+    %
+    % Returns:
+    %     dy: array containing the first derivative values.
 
     x_size = size(x, 2);
     y_size = size(y, 2);
@@ -24,7 +27,9 @@
     h = x(2) - x(1);
     n = x_size;
     dy = zeros(1, n);
+
     for i = 1:n
+
         if i == 1
             dy(i) = dy_end(h, y(i), y(i + 1), y(i + 2));
         elseif i == n
@@ -32,5 +37,7 @@
         else
             dy(i) = dy_mid(h, y(i - 1), y(i + 1));
         end
+
     end
+
 end
diff --git a/integration/composite2_simpson.m b/integration/composite2_simpson.m
index 4c41103..4d506b7 100644
--- a/integration/composite2_simpson.m
+++ b/integration/composite2_simpson.m
@@ -1,10 +1,12 @@
 function [xi] = composite2_simpson(x, y)
-% Calculate the integral from 1/3 Simpson's Rule
-% Inputs:
-%         x: Array containing x values
-%         y: Array containing y values
-% Outputs:
-%        xi: Integral value
+    % Calculate the integral from 1/3 Simpson's Rule.
+    %
+    % Args:
+    %     x: array containing x values.
+    %     y: array containing y values.
+    %
+    % Returns:
+    %     xi: integral value.
 
     x_size = size(x, 2);
     y_size = size(y, 2);
@@ -20,11 +22,13 @@
     sum_even = 0;
 
     for i = 2:(n - 1)
+
         if mod(i, 2) == 0
             sum_even = sum_even + y(i);
         else
             sum_odd = sum_odd + y(i);
         end
+
     end
 
     xi = h / 3 * (y(1) + 2 * sum_even + 4 * sum_odd + y(n));
diff --git a/integration/composite2_trapezoidal.m b/integration/composite2_trapezoidal.m
index 6748f9f..9d50b48 100644
--- a/integration/composite2_trapezoidal.m
+++ b/integration/composite2_trapezoidal.m
@@ -1,10 +1,12 @@
 function [xi] = composite2_trapezoidal(x, y)
-% Calculate the integral from Trapezoidal Rule
-% Inputs:
-%         x: Array containing x values
-%         y: Array containing y values
-% Outputs:
-%        xi: Integral value
+    % Calculate the integral from Trapezoidal Rule.
+    %
+    % Args:
+    %     x: array containing x values.
+    %     y: array containing y values.
+    %
+    % Returns:
+    %     xi: integral value.
 
     x_size = size(x, 2);
     y_size = size(y, 2);
diff --git a/integration/composite_simpson.m b/integration/composite_simpson.m
index 29adaf0..5b93f99 100644
--- a/integration/composite_simpson.m
+++ b/integration/composite_simpson.m
@@ -1,12 +1,14 @@
 function [xi] = composite_simpson(f, b, a, n)
-% Calculate the integral from 1/3 Simpson's Rule
-% Inputs:
-%         f: Function f(x)
-%         a: Initial point
-%         b: End point
-%         n: Number of intervals
-% Outputs:
-%        xi: Integral value
+    % Calculate the integral from 1/3 Simpson's Rule.
+    %
+    % Args:
+    %     f: function f(x).
+    %     a: initial point.
+    %     b: end point.
+    %     n: number of intervals.
+    %
+    % Returns:
+    %     xi: integral value.
 
     h = (b - a) / n;
 
@@ -15,11 +17,13 @@
 
     for i = 1:(n - 1)
         x = a + i * h;
+
         if mod(i, 2) == 0
             sum_even = sum_even + f(x);
         else
             sum_odd = sum_odd + f(x);
         end
+
     end
 
     xi = h / 3 * (f(a) + 2 * sum_even + 4 * sum_odd + f(b));
diff --git a/integration/composite_trapezoidal.m b/integration/composite_trapezoidal.m
index eaa7498..b3446d8 100644
--- a/integration/composite_trapezoidal.m
+++ b/integration/composite_trapezoidal.m
@@ -1,12 +1,14 @@
 function [xi] = composite_trapezoidal(f, b, a, n)
-% Calculate the integral from Trapezoidal Rule
-% Inputs:
-%         f: Function f(x)
-%         a: Initial point
-%         b: End point
-%         n: Number of intervals
-% Outputs:
-%        xi: Integral value
+    % Calculate the integral from Trapezoidal Rule.
+    %
+    % Args:
+    %     f: function f(x).
+    %     a: initial point.
+    %     b: end point.
+    %     n: number of intervals.
+    %
+    % Returns:
+    %     xi: integral value.
 
     h = (b - a) / n;
 
diff --git a/interpolation/lagrange.m b/interpolation/lagrange.m
index 601aae2..69ae71f 100644
--- a/interpolation/lagrange.m
+++ b/interpolation/lagrange.m
@@ -1,22 +1,30 @@
 function [y_int] = lagrange(x, y, x_int)
-% Interpolates a value using Lagrange polynomial
-% Inputs:
-%         x: Array containing x values
-%         y: Array containing y values
-%     x_int: Value to interpolate
-% Outputs:
-%     y_int: Interpolated value
+    % Interpolates a value using Lagrange polynomial.
+    %
+    % Args:
+    %     x: array containing x values.
+    %     y: array containing y values.
+    %     x_int: value to interpolate.
+    %
+    % Returns:
+    %     y_int: interpolated value.
 
     n = size(x, 2);
 
     y_int = 0;
+
     for i = 1:n
         p = y(i);
+
         for j = 1:n
+
             if i ~= j
                 p = p * ((x_int - x(j)) / (x(i) - x(j)));
             end
+
         end
+
         y_int = y_int + p;
     end
+
 end
diff --git a/interpolation/neville.m b/interpolation/neville.m
index 71db14c..ee1bee2 100644
--- a/interpolation/neville.m
+++ b/interpolation/neville.m
@@ -1,22 +1,26 @@
 function [y_int, q] = neville(x, y, x_int)
-% Interpolates a value using Neville polynomial
-% Inputs:
-%         x: Array containing x values
-%         y: Array containing y values
-%     x_int: Value to interpolate
-% Outputs:
-%     y_int: Interpolated value
-%         q: Coefficients matrix
+    % Interpolates a value using Neville polynomial.
+    %
+    % Args:
+    %     x: array containing x values.
+    %     y: array containing y values.
+    %     x_int: value to interpolate.
+    %
+    % Returns:
+    %     y_int: interpolated value.
+    %     q: coefficients matrix.
 
     n = size(x, 2);
     q = zeros(n, n - 1);
-    q = [y' q];  % Insert 'y' in the first column of the matrix 'q'
+    q = [y' q]; % Insert 'y' in the first column of the matrix 'q'
 
     for i = 2:n
+
         for j = 2:i
             q(i, j) = (x_int - x(i - j + 1)) * q(i, j - 1) - (x_int - x(i)) * q(i - 1, j - 1);
             q(i, j) = q(i, j) / (x(i) - x(i - j + 1));
         end
+
     end
 
     y_int = q(n, n);
diff --git a/linear_systems/backward_substitution.m b/linear_systems/backward_substitution.m
index 9cfa0c7..c1597b1 100644
--- a/linear_systems/backward_substitution.m
+++ b/linear_systems/backward_substitution.m
@@ -1,10 +1,12 @@
 function [x] = backward_substitution(u, d)
-% Solve the upper linear system ux=d
-% Inputs:
-%         u: Upper triangular matrix
-%         d: Array containing d values
-% Outputs:
-%         x: Solution of linear system
+    % Solve the upper linear system ux=d.
+    %
+    % Args:
+    %     upper: upper triangular matrix.
+    %     d: array containing d values.
+    %
+    % Returns:
+    %     x: solution of linear system.
 
     [n, m] = size(u);
 
@@ -15,6 +17,7 @@
     x = zeros(1, n);
 
     for i = n:-1:1
+
         if (u(i, i) == 0)
             error('Error: Matrix "u" is singular.')
         end
@@ -22,4 +25,5 @@
         x(i) = d(i) / u(i, i);
         d(1:i - 1) = d(1:i - 1) - u(1:i - 1, i) * x(i);
     end
+
 end
diff --git a/linear_systems/forward_substitution.m b/linear_systems/forward_substitution.m
index 3373dba..6105651 100644
--- a/linear_systems/forward_substitution.m
+++ b/linear_systems/forward_substitution.m
@@ -1,10 +1,12 @@
 function [x] = forward_substitution(l, c)
-% Solve the lower linear system lx=c
-% Inputs:
-%         l: Lower triangular matrix
-%         c: Array containing c values
-% Outputs:
-%         x: Solution of linear system
+    % Solve the lower linear system lx=c.
+    %
+    % Args:
+    %     lower: lower triangular matrix.
+    %     c: array containing c values.
+    %
+    % Returns:
+    %     x: solution of linear system.
 
     [n, m] = size(l);
 
@@ -15,6 +17,7 @@
     x = zeros(1, n);
 
     for i = 1:n
+
         if (l(i, i) == 0)
             error('Error: Matrix "l" is singular.')
         end
@@ -22,4 +25,5 @@
         x(i) = c(i) / l(i, i);
         c(i + 1:n) = c(i + 1:n) - l(i + 1:n, i) * x(i);
     end
+
 end
diff --git a/linear_systems/gauss_elimination_pp.m b/linear_systems/gauss_elimination_pp.m
index 8b3b875..48decef 100644
--- a/linear_systems/gauss_elimination_pp.m
+++ b/linear_systems/gauss_elimination_pp.m
@@ -1,10 +1,15 @@
 function [a] = gauss_elimination_pp(a, b)
-% Gaussian Elimination with Partial Pivoting - Calculate the upper upper triangular matrix from linear system Ax=b (do a row reduction)
-% Inputs:
-%         a: Matrix A from system Ax=b
-%         b: Array containing b values
-% Outputs:
-%         a: Augmented upper triangular matrix
+    % Gaussian Elimination with Partial Pivoting.
+    %
+    % Calculate the upper triangular matrix from linear system Ax=b (do a row
+    % reduction).
+    %
+    % Args:
+    %     a: matrix A from system Ax=b.
+    %     b: array containing b values.
+    %
+    % Returns:
+    %     a: augmented upper triangular matrix.
 
     [n, m] = size(a);
 
@@ -20,10 +25,12 @@
 
         % Comparison to select the pivot
         for j = (i + 1):n
+
             if abs(a(j, i)) > abs(a(i, i))
                 % Swap rows
                 a([i j], :) = a([j i], :);
             end
+
         end
 
         % Cheking for nullity of the pivots
@@ -34,19 +41,23 @@
         if p == n + 1
             warning('Info: No unique solution.');
         else
+
             if p ~= i
                 % Swap rows
                 a([i p], :) = a([p i], :);
             end
+
         end
 
         for j = (i + 1):n
             a(j, :) = a(j, :) - a(i, :) * (a(j, i) / a(i, i));
         end
+
     end
 
     % Checking for nonzero of last entry
     if a(n, n) == 0
         warning('Info: No unique solution.');
     end
+
 end
diff --git a/linear_systems_iterative/gauss_seidel.m b/linear_systems_iterative/gauss_seidel.m
index 4cacec6..38ad44a 100644
--- a/linear_systems_iterative/gauss_seidel.m
+++ b/linear_systems_iterative/gauss_seidel.m
@@ -1,14 +1,16 @@
 function [x, iter] = gauss_seidel(a, b, x0, tol, iter_max)
-% Gauss-Seidel method: solve Ax = b given an initial approximation x0
-% Inputs:
-%         a: Matrix A from system Ax=b
-%         b: Array containing b values
-%        x0: Initial approximation of solution
-%       tol: Tolerance
-%  iter_max: Maximum number of iterations
-% Outpus:
-%         x: Solution of linear system
-%      iter: Used iterations
+    % Gauss-Seidel method: solve Ax = b given an initial approximation x0.
+    %
+    % Args:
+    %     a: matrix A from system Ax=b.
+    %     b: array containing b values.
+    %     x0: initial approximation of solution.
+    %     tol: tolerance.
+    %     iter_max: maximum number of iterations.
+    %
+    % Returns:
+    %     x: solution of linear system.
+    %     iter: used iterations.
 
     % L and U matrices
     l = tril(a);
@@ -16,11 +18,13 @@
 
     % Iterative process
     for iter = 1:iter_max
-        x = l \ (b - u * x0);  % "A\B" is the same as "INV(A)*B"
+        x = l \ (b - u * x0); % "A\B" is the same as "INV(A)*B"
 
         if norm(x - x0, inf) / norm(x, inf) <= tol
             break;
         end
+
         x0 = x;
     end
+
 end
diff --git a/linear_systems_iterative/jacobi.m b/linear_systems_iterative/jacobi.m
index 039c5cf..0bc0388 100644
--- a/linear_systems_iterative/jacobi.m
+++ b/linear_systems_iterative/jacobi.m
@@ -1,14 +1,16 @@
 function [x, iter] = jacobi(a, b, x0, tol, iter_max)
-% Jacobi method: solve Ax = b given an initial approximation x0
-% Inputs:
-%         a: Matrix A from system Ax=b
-%         b: Array containing b values
-%        x0: Initial approximation of solution
-%       tol: Tolerance
-%  iter_max: Maximum number of iterations
-% Outpus:
-%         x: Solution of linear system
-%      iter: Used iterations
+    % Jacobi method: solve Ax = b given an initial approximation x0.
+    %
+    % Args:
+    %     a: matrix A from system Ax=b.
+    %     b: array containing b values.
+    %     x0: initial approximation of solution.
+    %     tol: tolerance.
+    %     iter_max: maximum number of iterations.
+    %
+    % Returns:
+    %     x: solution of linear system.
+    %     iter: used iterations.
 
     % D and M matrices
     d = diag(diag(a));
@@ -16,11 +18,13 @@
 
     % Iterative process
     for iter = 1:iter_max
-        x = d \ (b - m * x0);  % "A\B" is the same as "INV(A)*B"
+        x = d \ (b - m * x0); % "A\B" is the same as "INV(A)*B"
 
         if norm(x - x0, inf) / norm(x, inf) <= tol
             break;
         end
+
         x0 = x;
     end
+
 end
diff --git a/main.m b/main.m
index 5e894fa..2b8f44e 100644
--- a/main.m
+++ b/main.m
@@ -1,7 +1,6 @@
 %% Numerical methods implementation in MATLAB.
 %
 % Author: Cristiano Nunes <cfgnunes@gmail.com>
-% Date: 2017-10-29
 
 clc;
 close all;
@@ -18,152 +17,157 @@
 addpath solutions
 
 % Bisection method (find roots of an equation)
-%     Pros:
-%         It is a reliable method with guaranteed convergence;
-%         It is a simple method that does the search of the root by means of a binary search;
-%         There is no need to calculate the derivative of the function.
-%     Cons:
-%         Slow convergence;
-%         It is necessary to enter a search interval [a, b];
-%         The interval reported must have a signal exchange, f (a) * f (b) <0.
-disp('> Running Solutions: Bisection method');
-f = @(x) (4 * x ^ 3 + x + cos(x) - 10);
-tol = 10 ^ -5;
+%   Pros:
+%       It is a reliable method with guaranteed convergence;
+%       It is a simple method that searches for the root employing a
+%           binary search;
+%       There is no need to calculate the derivative of the function.
+%   Cons:
+%       Slow convergence;
+%       It is necessary to enter a search interval [a, b];
+%       The interval reported must have a signal exchange, f (a) * f (b)<0.
+disp('> Run an example "Solutions: Bisection method".')
+f = @(x) (4 * x^3 + x + cos(x) - 10);
+tol = 10^ - 5;
 iter_max = 100;
 a = 1.0;
 b = 2.0;
 [root, iter, converged] = bisection(f, a, b, tol, iter_max)
 
 % Newton method (find roots of an equation)
-%     Pros:
-%         It is a fast method.
-%     Cons:
-%         It may diverge;
-%         It is necessary to calculate the derivative of the function;
-%         It is necessary to give an initial x0 value where f'(x0) must be nonzero.
-disp('> Running Solutions: Newton method');
-f = @(x) (4 * x ^ 3 + x + cos(x) - 10);
-df = @(x) (12 * x ^ 2 + 1 - sin(x));
-tol = 10 ^ -5;
+%   Pros:
+%       It is a fast method.
+%    Cons:
+%       It may diverge;
+%       It is necessary to calculate the derivative of the function;
+%       It is necessary to give an initial x0 value where
+%           f'(x0) must be nonzero.
+disp('> Run an example "Solutions: Newton method".')
+f = @(x) (4 * x^3 + x + cos(x) - 10);
+df = @(x) (12 * x^2 + 1 - sin(x));
+tol = 10^ - 5;
 iter_max = 100;
 x0 = 2.0;
 [root, iter, converged] = newton(f, df, x0, tol, iter_max)
 
 % Secant method (find roots of an equation)
-%     Pros:
-%         It is a fast method (slower than Newton's method);
-%         It is based on the Newton method, but does not need the derivative of the function.
-%     Cons:
-%         It may diverge if the function is not approximately linear in the range containing the root;
-%         It is necessary to give two points 'a' and 'b' where f(a)-f(b) must be nonzero.
-disp('> Running Solutions: Secant method');
-f = @(x) (4 * x ^ 3 + x + cos(x) - 10);
-tol = 10 ^ -5;
+%   Pros:
+%       It is a fast method (slower than Newton's method);
+%       It is based on the Newton method but does not need the derivative
+%           of the function.
+%   Cons:
+%       It may diverge if the function is not approximately linear in the
+%           range containing the root;
+%       It is necessary to give two points 'a' and 'b' where
+%           f(a)-f(b) must be nonzero.
+disp('> Run an example "Solutions: Secant method".')
+f = @(x) (4 * x^3 + x + cos(x) - 10);
+tol = 10^ - 5;
 iter_max = 100;
 a = 1.0;
 b = 2.0;
 [root, iter, converged] = secant(f, a, b, tol, iter_max)
 
-disp('> Running Interpolation: Lagrange method');
-x = [2 11 / 4 4];
-y = [1 / 2 4 / 11 1 / 4];
+disp('> Run an example "Interpolation: Lagrange method".')
+x = [2 11/4 4];
+y = [1/2 4/11 1/4];
 x_int = 3;
 [y_int] = lagrange(x, y, x_int)
 
-disp('> Running Interpolation: Neville method');
+disp('> Run an example "Interpolation: Neville method".')
 x = [1.0 1.3 1.6 1.9 2.2];
 y = [0.7651977 0.6200860 0.4554022 0.2818186 0.1103623];
 x_int = 1.5;
 [y_int, q] = neville(x, y, x_int)
 
-disp('> Running Polynomials: Briot-Ruffini method');
+disp('> Run an example "Polynomials: Briot-Ruffini method".')
 root = -2;
 a = [2 0 -3 3 -4];
 [b, rest] = briot_ruffini(root, a)
 
-disp('> Running Polynomials: Newtons Divided-Difference method');
+disp('> Run an example "Polynomials: Newtons Divided-Difference method".')
 x = [1.0 1.3 1.6 1.9 2.2];
 y = [0.7651977 0.6200860 0.4554022 0.2818186 0.1103623];
 [f] = newton_divided_difference(x, y)
 
-disp('> Running Differentiation: Backward-difference method')
+disp('> Run an example "Differentiation: Backward-difference method".')
 x = [0.0 0.2 0.4];
 y = [0.00000 0.74140 1.3718];
 [dy] = derivative_backward_difference(x, y)
 
-disp('> Running Differentiation: Three-Point method')
+disp('> Run an example "Differentiation: Three-Point method".')
 x = [1.1 1.2 1.3 1.4];
 y = [9.025013 11.02318 13.46374 16.44465];
 [dy] = derivative_three_point(x, y)
 
-disp('> Running Differentiation: Five-Point method')
+disp('> Run an example "Differentiation: Five-Point method".')
 x = [2.1 2.2 2.3 2.4 2.5 2.6];
 y = [-1.709847 -1.373823 -1.119214 -0.9160143 -0.7470223 -0.6015966];
 [dy] = derivative_five_point(x, y)
 
-disp('> Running Integration: Trapezoidal Rule')
+disp('> Run an example "Integration: Trapezoidal Rule".')
 x = [0, 6, 12, 18, 24, 30, 36, 42, 48, 54, 60, 66, 72, 78, 84];
 y = [124, 134, 148, 156, 147, 133, 121, 109, 99, 85, 78, 89, 104, 116, 123];
 [xi] = composite2_trapezoidal(x, y)
 
-disp('> Running Integration: Trapezoidal Rule')
-f = @(x) (x ^ 2 * log(x ^ 2 + 1));
+disp('> Run an example "Integration: Trapezoidal Rule".')
+f = @(x) (x^2 * log(x^2 + 1));
 a = 0.0;
 b = 2.0;
 h = 0.25;
 n = (b - a) / h;
 [xi] = composite_trapezoidal(f, b, a, n)
 
-disp('> Running Integration: Composite 1/3 Simpsons Rule')
+disp('> Run an example "Integration: Composite 1/3 Simpsons Rule".')
 x = [0, 6, 12, 18, 24, 30, 36, 42, 48, 54, 60, 66, 72, 78, 84];
 y = [124, 134, 148, 156, 147, 133, 121, 109, 99, 85, 78, 89, 104, 116, 123];
 [xi] = composite2_simpson(x, y)
 
-disp('> Running Integration: Composite 1/3 Simpsons Rule')
-f = @(x) (x ^ 2 * log(x ^ 2 + 1));
+disp('> Run an example "Integration: Composite 1/3 Simpsons Rule".')
+f = @(x) (x^2 * log(x^2 + 1));
 a = 0.0;
 b = 2.0;
 h = 0.25;
 n = (b - a) / h;
 [xi] = composite_simpson(f, b, a, n)
 
-disp('> Running ODE: Euler method')
-f = @(x, y) (y - x ^ 2 + 1);
+disp('> Run an example "ODE: Euler method".')
+f = @(x, y) (y - x^2 + 1);
 a = 0.0;
 b = 2.0;
 n = 10;
 ya = 0.5;
 [vx, vy] = euler(f, a, b, n, ya);
 
-disp('> Running ODE: Taylor (Order Two) method')
-f = @(x, y) (y - x ^ 2 + 1);
-df1 = @(x, y) (y - x ^ 2 + 1 - 2 * x);
+disp('> Run an example "ODE: Taylor (Order Two) method".')
+f = @(x, y) (y - x^2 + 1);
+df1 = @(x, y) (y - x^2 + 1 - 2 * x);
 a = 0.0;
 b = 2.0;
 n = 10;
 ya = 0.5;
 [vx, vy] = taylor2(f, df1, a, b, n, ya);
 
-disp('> Running ODE: Taylor (Order Four) method')
-f = @(x, y) (y - x ^ 2 + 1);
-df1 = @(x, y) (y - x ^ 2 + 1 - 2 * x);
-df2 = @(x, y) (y - x ^ 2 + 1 - 2 * x - 2);
-df3 = @(x, y) (y - x ^ 2 + 1 - 2 * x - 2);
+disp('> Run an example "ODE: Taylor (Order Four) method".')
+f = @(x, y) (y - x^2 + 1);
+df1 = @(x, y) (y - x^2 + 1 - 2 * x);
+df2 = @(x, y) (y - x^2 + 1 - 2 * x - 2);
+df3 = @(x, y) (y - x^2 + 1 - 2 * x - 2);
 a = 0.0;
 b = 2.0;
 n = 10;
 ya = 0.5;
 [vx, vy] = taylor4(f, df1, df2, df3, a, b, n, ya);
 
-disp('> Running ODE: Runge-Kutta (Order Four) method')
-f = @(x, y) (y - x ^ 2 + 1);
+disp('> Run an example "ODE: Runge-Kutta (Order Four) method".')
+f = @(x, y) (y - x^2 + 1);
 a = 0.0;
 b = 2.0;
 n = 10;
 ya = 0.5;
 [vx, vy] = rk4(f, a, b, n, ya);
 
-disp('> Running ODE: Runge-Kutta (Order Four) method for systems of differential equations')
+disp('> Run an example "ODE: Runge-Kutta (Order Four) method for systems of differential equations".')
 m = 2;
 f = cell(m, 1);
 f{1} = @(x, y) (- 4 * y(1) + 3 * y(2) + 6);
@@ -177,33 +181,33 @@
 ya(2) = 0.0;
 [vx, vy] = rk4_system(f, a, b, n, ya)
 
-disp('> Running Linear Systems: Gaussian Elimination')
+disp('> Run an example "Linear Systems: Gaussian Elimination".')
 a = [1 -1 2 -1; 2 -2 3 -3; 1 1 1 0; 1 -1 4 3];
 b = [-8 -20 -2 4];
 [a] = gauss_elimination_pp(a, b')
 
-disp('> Running Linear Systems: Backward Substitution')
-u = a(:,1:end-1);
-d = a(:,end);
+disp('> Run an example "Linear Systems: Backward Substitution".')
+u = a(:, 1:end - 1);
+d = a(:, end);
 [x] = backward_substitution(u, d)
 
-disp('> Running Linear Systems: Forward Substitution')
+disp('> Run an example "Linear Systems: Forward Substitution".')
 l = [3 0 0 0; -1 1 0 0; 3 -2 -1 0; 1 -2 6 2];
 c = [5 6 4 2];
 [a] = forward_substitution(l, c')
 
-disp('> Running Iteractive Linear Systems: Jacobi')
+disp('> Run an example "Iteractive Linear Systems: Jacobi".')
 a = [10 -1 2 0; -1 11 -1 3; 2 -1 10 -1; 0 3 -1 8];
 b = [6 25 -11 15];
 x0 = [0 0 0 0];
-tol = 10 ^ -3;
+tol = 10^ - 3;
 iter_max = 10;
 [x, iter] = jacobi(a, b', x0', tol, iter_max)
 
-disp('> Running Iteractive Linear Systems: Gauss-Seidel')
+disp('> Run an example "Iteractive Linear Systems: Gauss-Seidel".')
 a = [10 -1 2 0; -1 11 -1 3; 2 -1 10 -1; 0 3 -1 8];
 b = [6 25 -11 15];
 x0 = [0 0 0 0];
-tol = 10 ^ -3;
+tol = 10^ - 3;
 iter_max = 10;
 [x, iter] = gauss_seidel(a, b', x0', tol, iter_max)
diff --git a/ode/euler.m b/ode/euler.m
index 15acbb3..8b25ceb 100644
--- a/ode/euler.m
+++ b/ode/euler.m
@@ -1,14 +1,18 @@
 function [vx, vy] = euler(f, a, b, n, ya)
-% Calculate the solution of the initial-value problem from Euler method
-% Inputs:
-%         f: Function f(x)
-%         a: Initial point
-%         b: End point
-%         n: Number of intervals
-%        ya: Initial value
-% Outputs:
-%       vx: Array containing x values
-%       vy: Array containing y values (solution of IVP)
+    % Calculate the solution of the initial-value problem (IVP).
+    %
+    % Solve the IVP from Euler method.
+    %
+    % Args:
+    %     f: function f(x).
+    %     a: initial point.
+    %     b: end point.
+    %     n: number of intervals.
+    %     ya: initial value.
+    %
+    % Returns:
+    %     vx: array containing x values.
+    %     vy: array containing y values (solution of IVP).
 
     vx = zeros(1, n + 1);
     vy = zeros(1, n + 1);
@@ -32,4 +36,5 @@
         vx(i + 1) = x;
         vy(i + 1) = y;
     end
+
 end
diff --git a/ode/rk4.m b/ode/rk4.m
index 8a1f3c4..47417ca 100644
--- a/ode/rk4.m
+++ b/ode/rk4.m
@@ -1,14 +1,18 @@
 function [vx, vy] = rk4(f, a, b, n, ya)
-% Calculate the solution of the initial-value problem from Runge-Kutta (Order Four) method
-% Inputs:
-%         f: Function f(x)
-%         a: Initial point
-%         b: End point
-%         n: Number of intervals
-%        ya: Initial value
-% Outputs:
-%       vx: Array containing x values
-%       vy: Array containing y values (solution of IVP)
+    % Calculate the solution of the initial-value problem (IVP).
+    %
+    % Solve the IVP from Runge-Kutta (Order Four) method.
+    %
+    % Args:
+    %     f: function f(x).
+    %     a: initial point.
+    %     b: end point.
+    %     n: number of intervals.
+    %     ya: initial value.
+    %
+    % Returns:
+    %     vx: array containing x values.
+    %     vy: array containing y values (solution of IVP).
 
     vx = zeros(1, n + 1);
     vy = zeros(1, n + 1);
@@ -35,4 +39,5 @@
         vx(i + 1) = x;
         vy(i + 1) = y;
     end
+
 end
diff --git a/ode/rk4_system.m b/ode/rk4_system.m
index eb116eb..ea29598 100644
--- a/ode/rk4_system.m
+++ b/ode/rk4_system.m
@@ -1,14 +1,18 @@
 function [vx, vy] = rk4_system(f, a, b, n, ya)
-% Calculate the solution of systems of differential equations from Runge-Kutta (Order Four) method
-% Inputs:
-%         f: Array of functions f(x)
-%         a: Initial point
-%         b: End point
-%         n: Number of intervals
-%        ya: Array of initial values
-% Outputs:
-%       vx: Array containing x values
-%       vy: Array containing y values (solution of IVP)
+    % Calculate the solution of systems of differential equations.
+    %
+    % Solve from Runge-Kutta (Order Four) method.
+    %
+    % Args:
+    %     f: array of functions f(x).
+    %     a: initial point.
+    %     b: end point.
+    %     n: number of intervals.
+    %     ya: array of initial values.
+    %
+    % Returns:
+    %     vx: array containing x values.
+    %     vy: array containing y values (solution of IVP).
 
     m = size(f, 1);
 
@@ -52,4 +56,5 @@
         vx(i + 1) = x;
         vy(:, i + 1) = y;
     end
+
 end
diff --git a/ode/taylor2.m b/ode/taylor2.m
index 2781918..b529665 100644
--- a/ode/taylor2.m
+++ b/ode/taylor2.m
@@ -1,15 +1,19 @@
 function [vx, vy] = taylor2(f, df1, a, b, n, ya)
-% Calculate the solution of the initial-value problem from Taylor (Order Two) method
-% Inputs:
-%         f: Function f(x)
-%       df1: 1's derivative of function f(x)
-%         a: Initial point
-%         b: End point
-%         n: Number of intervals
-%        ya: Initial value
-% Outputs:
-%       vx: Array containing x values
-%       vy: Array containing y values (solution of IVP)
+    % Calculate the solution of the initial-value problem (IVP).
+    %
+    % Solve the IVP from Taylor (Order Two) method.
+    %
+    % Args:
+    %     f: function f(x).
+    %     df1: 1's derivative of function f(x).
+    %     a: initial point.
+    %     b: end point.
+    %     n: number of intervals.
+    %     ya: initial value.
+    %
+    % Returns:
+    %     vx: array containing x values.
+    %     vy: array containing y values (solution of IVP).
 
     vx = zeros(1, n + 1);
     vy = zeros(1, n + 1);
@@ -31,4 +35,5 @@
         vx(i + 1) = x;
         vy(i + 1) = y;
     end
+
 end
diff --git a/ode/taylor4.m b/ode/taylor4.m
index 69fb553..8ea87e8 100644
--- a/ode/taylor4.m
+++ b/ode/taylor4.m
@@ -1,17 +1,21 @@
 function [vx, vy] = taylor4(f, df1, df2, df3, a, b, n, ya)
-% Calculate the solution of the initial-value problem from Taylor (Order Four) method
-% Inputs:
-%         f: Function f(x)
-%       df1: 1's derivative of function f(x)
-%       df2: 2's derivative of function f(x)
-%       df3: 3's derivative of function f(x)
-%         a: Initial point
-%         b: End point
-%         n: Number of intervals
-%        ya: Initial value
-% Outputs:
-%       vx: Array containing x values
-%       vy: Array containing y values (solution of IVP)
+    % Calculate the solution of the initial-value problem (IVP).
+    %
+    % Solve the IVP from Taylor (Order Four) method.
+    %
+    % Args:
+    %     f: function f(x).
+    %     df1: 1's derivative of function f(x).
+    %     df2: 2's derivative of function f(x).
+    %     df3: 3's derivative of function f(x).
+    %     a: initial point.
+    %     b: end point.
+    %     n: number of intervals.
+    %     ya: initial value.
+    %
+    % Returns:
+    %     vx: array containing x values.
+    %     vy: array containing y values (solution of IVP).
 
     vx = zeros(1, n + 1);
     vy = zeros(1, n + 1);
@@ -26,11 +30,12 @@
     fprintf('i: %.3d\t x:%.4f\t y:%.4f\t\n', 0, x, y);
 
     for i = 1:n
-        y = y + h * (f(x, y) + 0.5 * h * df1(x, y) + (h ^ 2 / 6) * df2(x, y) + (h ^ 3 / 24) * df3(x, y));
+        y = y + h * (f(x, y) + 0.5 * h * df1(x, y) + (h^2/6) * df2(x, y) + (h^3/24) * df3(x, y));
         x = a + i * h;
 
         fprintf('i: %.3d\t x:%.4f\t y:%.4f\t\n', i, x, y);
         vx(i + 1) = x;
         vy(i + 1) = y;
     end
+
 end
diff --git a/polynomials/briot_ruffini.m b/polynomials/briot_ruffini.m
index af17685..19fc39a 100644
--- a/polynomials/briot_ruffini.m
+++ b/polynomials/briot_ruffini.m
@@ -1,12 +1,15 @@
 function [b, rest] = briot_ruffini(root, a)
-% Divides a polynomial by another polynomial in the format (x-root)
-% P(x) = Q(x) * (x-root) + rest
-% Inputs:
-%         a: Array containing the coefficients of the input polynomial
-%      root: One of the polynomial roots
-% Outpus:
-%         b: Array containing the coefficients of the output polynomial
-%      rest: Polynomial division Rest
+    % Divides a polynomial by another polynomial.
+    %
+    % The format is: P(x) = Q(x) * (x-root) + rest.
+    %
+    % Args:
+    %     a: array containing the coefficients of the input polynomial.
+    %     root: one of the polynomial roots.
+    %
+    % Returns:
+    %     b: array containing the coefficients of the output polynomial.
+    %     rest: polynomial division Rest.
 
     n = size(a, 2) - 1;
     b = zeros(1, n);
diff --git a/polynomials/newton_divided_difference.m b/polynomials/newton_divided_difference.m
index e98cd10..d0a30d3 100644
--- a/polynomials/newton_divided_difference.m
+++ b/polynomials/newton_divided_difference.m
@@ -1,24 +1,31 @@
 function [f] = newton_divided_difference(x, y)
-% Find the coefficients of Newton's divided difference and the Newton's polynomial
-% Inputs:
-%         x: Array containing x values
-%         y: Array containing y values
-% Outpus:
-%         f: Array containing Newton's divided difference coefficients
+    % Find the coefficients of Newton's divided difference.
+    %
+    % Also findthe Newton's polynomial.
+    %
+    % Args:
+    %     x: array containing x values.
+    %     y: array containing y values.
+    %
+    % Returns:
+    %     f: array containing Newton's divided difference coefficients.
 
     n = size(x, 2);
     q = zeros(n, n - 1);
-    q = [y' q];  % Insert 'y' in the first column of the matrix 'q'
+    q = [y' q]; % Insert 'y' in the first column of the matrix 'q'
 
     for i = 2:n
+
         for j = 2:i
             q(i, j) = q(i, j - 1) - q(i - 1, j - 1);
             q(i, j) = q(i, j) / (x(i) - x(i - j + 1));
         end
+
     end
 
     % Copy the diagonal values of the matrix q to the vector f
     f = zeros(1, n);
+
     for i = 1:n
         f(i) = q(i, i);
     end
@@ -26,12 +33,16 @@
     % Prints the polynomial
     disp('The polynomial is:')
     fprintf('P(x)=%+.4f', f(1));
+
     for i = 2:n
         fprintf('%+.4f', f(i));
+
         for j = 2:i
-            fprintf('(x%+.4f)', x(j - 1) * - 1);
+            fprintf('(x%+.4f)', x(j - 1) *- 1);
         end
+
     end
+
     fprintf('\n');
 
 end
diff --git a/solutions/bisection.m b/solutions/bisection.m
index 9c77ace..6f02667 100644
--- a/solutions/bisection.m
+++ b/solutions/bisection.m
@@ -1,15 +1,17 @@
 function [root, iter, converged] = bisection(f, a, b, tol, iter_max)
-% Calculates the root of an equation by Bisection method
-% Inputs:
-%         f: Function f(x)
-%         a: Lower limit
-%         b: Upper limit
-%       tol: Tolerance
-%  iter_max: Maximum number of iterations
-% Outpus:
-%      root: Root value
-%      iter: Used iterations
-% converged: Found the root
+    % Calculate the root of an equation by Bisection method.
+    %
+    % Args:
+    %     f: function f(x).
+    %     a: lower limit.
+    %     b: upper limit.
+    %     tol: tolerance.
+    %     iter_max: maximum number of iterations.
+    %
+    % Returns:
+    %     root: root value.
+    %     iter: used iterations.
+    %     converged: found the root.
 
     fa = f(a);
     fb = f(b);
@@ -48,4 +50,5 @@
         warning('Warning: The method did not converge.');
         converged = 0;
     end
+
 end
diff --git a/solutions/newton.m b/solutions/newton.m
index 4332a97..0474899 100644
--- a/solutions/newton.m
+++ b/solutions/newton.m
@@ -1,15 +1,17 @@
 function [root, iter, converged] = newton(f, df, x0, tol, iter_max)
-% Calculates the root of an equation by Newton method
-% Inputs:
-%         f: Function f(x)
-%        df: Derivative of function f(x)
-%        x0: Initial guess
-%       tol: Tolerance
-%  iter_max: Maximum number of iterations
-% Outpus:
-%      root: Root value
-%      iter: Used iterations
-% converged: Found the root
+    % Calculate the root of an equation by Newton method.
+    %
+    % Args:
+    %     f: function f(x).
+    %     df: derivative of function f(x).
+    %     x0: initial guess.
+    %     tol: tolerance.
+    %     iter_max: maximum number of iterations.
+    %
+    % Returns:
+    %     root: root value.
+    %     iter: used iterations.
+    %     converged: found the root.
 
     x = x0;
     fx = f(x);
@@ -19,7 +21,7 @@
     fprintf('iter: %.3d\t x: %.4f\t dfx: %.4f\t fx: %.4f\n', iter, x, dfx, fx);
 
     for iter = 1:iter_max
-        deltaX = - fx / dfx;
+        deltaX = -fx / dfx;
         x = x + deltaX;
         fx = f(x);
         dfx = df(x);
@@ -29,6 +31,7 @@
         if (abs(deltaX) <= tol && abs(fx) <= tol) || dfx == 0
             break;
         end
+
     end
 
     root = x;
@@ -39,4 +42,5 @@
         warning('Warning: The method did not converge.');
         converged = 0;
     end
+
 end
diff --git a/solutions/secant.m b/solutions/secant.m
index 99236ef..1c3ec11 100644
--- a/solutions/secant.m
+++ b/solutions/secant.m
@@ -1,15 +1,17 @@
 function [root, iter, converged] = secant(f, a, b, tol, iter_max)
-% Calculates the root of an equation by Secant method
-% Inputs:
-%         f: Function f(x)
-%         a: Lower limit
-%         b: Upper limit
-%       tol: Tolerance
-%  iter_max: Maximum number of iterations
-% Outpus:
-%      root: Root value
-%      iter: Used iterations
-% converged: Found the root
+    % Calculate the root of an equation by Secant method.
+    %
+    % Args:
+    %     f: function f(x).
+    %     a: lower limit.
+    %     b: upper limit.
+    %     tol: tolerance.
+    %     iter_max: maximum number of iterations.
+    %
+    % Returns:
+    %     root: root value.
+    %     iter: used iterations.
+    %     converged: found the root.
 
     fa = f(a);
     fb = f(b);
@@ -33,8 +35,9 @@
 
     x = b;
     fx = fb;
+
     for iter = 0:iter_max
-        deltaX = - fx / (fb - fa) * (b - a);
+        deltaX = -fx / (fb - fa) * (b - a);
         x = x + deltaX;
         fx = f(x);
         fprintf('iter: %.3d\t a: %.4f\t fa: %.4f\t b: %.4f\t fb: %.4f\t x: %.4f\t fx: %.4f\t deltaX: %.4f\n', iter, a, fa, b, fb, x, fx, deltaX);
@@ -48,6 +51,7 @@
         b = x;
         fb = fx;
     end
+
     root = x;
 
     if abs(deltaX) <= tol && abs(fx) <= tol
@@ -56,4 +60,5 @@
         warning('Warning: The method did not converge.');
         converged = 0;
     end
+
 end