diff --git a/non_linear_algebra/doc.txt b/non_linear_algebra/doc.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f951b578fe82e29638a9ca07fea1180fd7dadac8
--- /dev/null
+++ b/non_linear_algebra/doc.txt
@@ -0,0 +1,22 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser
+ * General Public License v2.1. See the file LICENSE in the top level
+ * directory for more details.
+ */
+
+/**
+ * @defgroup non_linear_algebra NON_LINEAR_ALGEBRA
+ *
+ * @brief Non-linear algebra operations
+ *
+ * @details The non-linear algebra module contains functions to solve multi-variant nonlinear
+ * equations as wells algorithms solving problems of regression smoothing and curve
+ * fitting. This module also enables enables the optimization of an approximate solution
+ * by using Non-linear Least Squares (NLS) methods such as modified Gauss--Newton (GN) or
+ * the Levenberg--Marquardt (LVM) algorithms.
+ *
+ * @author Zakaria Kasmi
+ */
diff --git a/non_linear_algebra/optimization/doc.txt b/non_linear_algebra/optimization/doc.txt
new file mode 100644
index 0000000000000000000000000000000000000000..6586aab419a7339ea355847f69cb4507a60309c8
--- /dev/null
+++ b/non_linear_algebra/optimization/doc.txt
@@ -0,0 +1,22 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser
+ * General Public License v2.1. See the file LICENSE in the top level
+ * directory for more details.
+ */
+
+/**
+ * @defgroup optimization OPTIMIZATION
+ * @ingroup non_linear_algebra
+ *
+ * @brief Solving problems of regression smoothing and curve fitting.
+ *
+ * @details This module also enables the optimization of an
+ * approximate solution by using Non-linear Least Squares
+ * (NLS) methods such as modified Gauss--Newton (GN) or
+ * the Levenberg--Marquardt (LVM) algorithms.
+ *
+ * @author Zakaria Kasmi
+ */
diff --git a/non_linear_algebra/optimization/include/levenberg_marquardt.h b/non_linear_algebra/optimization/include/levenberg_marquardt.h
new file mode 100644
index 0000000000000000000000000000000000000000..7d162a9a65a46b72f9717eaf9c8217a2ec7d38a9
--- /dev/null
+++ b/non_linear_algebra/optimization/include/levenberg_marquardt.h
@@ -0,0 +1,81 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup optimization
+ * @{
+ *
+ * @file
+ * @brief Implement the Levenberg--Marquardt (LVM) algorithm.
+ * @note This function is generally implemented.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * @author Naouar Guerchali
+ *
+ * @}
+ */
+
+#ifndef LEVENBERG_MARQUARDT_H_
+#define LEVENBERG_MARQUARDT_H_
+
+#include <inttypes.h>
+
+#include "matrix.h"
+#include "vector.h"
+
+/**
+ * @brief Implements the Levenberg--Marquardt (LVM) algorithm.
+ * @details The user should provide pointers to the error and Jacobian functions.
+ * @note This function is generally implemented.
+ *
+ * @param[in] f_length length of the error functions vector.
+ * @param[in] n length of the start vector.
+ * @param[in] x0_vec[] start vector.
+ * @param[in] data_vec[] data vector.
+ * @param[in] eps accuracy bound.
+ * @param[in] tau \f$ \tau \f$ factor.
+ * @param[in] beta0 \f$ \beta_0 \f$ factor.
+ * @param[in] beta1 \f$ \beta_1 \f$ factor.
+ * @param[in] max_iter_num maximal iteration number of the LVM algorithm.
+ * @param[out] est_x_vec[] estimated (optimized) vector.
+ * @param[in] (*get_f_error) pointer to the error function.
+ * @param[in] (*get_jacobian) pointer to the Jacobian matrix.
+ *
+ * @return required iteration number.
+ *
+ */
+uint8_t opt_levenberg_marquardt(uint8_t f_length, uint8_t n,
+ vector_t x0_vec[n],
+ vector_t data_vec[f_length],
+ matrix_t eps, matrix_t tau, matrix_t beta0,
+ matrix_t beta1,
+ uint8_t max_iter_num,
+ vector_t est_x_vec[n],
+ void (*get_f_error)(vector_t x0_vec[],
+ vector_t data_vec[],
+ vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x0_vec[],
+ matrix_t J[][n])
+ );
+
+/**
+ * @brief Compute the initial value \f$ \mu_0 \f$ of the Levenberg--Marquardt (LVM) algorithm.
+ * @details The user should provide a pointer to the matrix
+ * \f$ J_f^{T} J_{f} \f$.
+ *
+ * @param[in] n column number of the matrix \f$J_f^T J_f\f$.
+ * @param[in] tau \f$ \tau \f$ factor.
+ * @param[in] JTJ[][] pointer to the matrix \f$J_f^T J_f\f$.
+ *
+ * @return parameter \f$ \rho_{\mu} \f$
+ */
+matrix_t opt_levenberg_marquardt_get_mu0(uint8_t n, matrix_t tau,
+ matrix_t JTJ[][n]);
+
+#endif /* LEVENBERG_MARQUARDT_H_ */
diff --git a/non_linear_algebra/optimization/include/modified_gauss_newton.h b/non_linear_algebra/optimization/include/modified_gauss_newton.h
new file mode 100644
index 0000000000000000000000000000000000000000..4182f9128efb6832982561a03ef84ca74b837185
--- /dev/null
+++ b/non_linear_algebra/optimization/include/modified_gauss_newton.h
@@ -0,0 +1,63 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup optimization
+ * @{
+ *
+ * @file
+ * @brief Implement the Gauss--Newton algorithm.
+ * @note This function is generally implemented.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * @author Abdelmoumen Norrdine <a.norrdine@googlemail.com>
+ *
+ * @}
+ */
+
+#ifndef MODIFIED_GAUSS_NEWTON_H_
+#define MODIFIED_GAUSS_NEWTON_H_
+
+#include <inttypes.h>
+
+#include "matrix.h"
+#include "vector.h"
+
+/**
+ * @brief Implements the modified Gauss--Newton algorithm.
+ * @details The user should provide pointers to the error and Jacobian functions.
+ * @note This function is generally implemented.
+ *
+ * @param[in] f_length length of the error functions vector.
+ * @param[in] n length of the start vector.
+ * @param[in] x0_vec[] start vector.
+ * @param[in] data_vec[] data vector.
+ * @param[in] eps accuracy bound.
+ * @param[in] fmin termination tolerance on the error function.
+ * @param[in] max_iter_num maximal iteration number of the Gauss--Newton algorithm.
+ * @param[out] est_x_vec[] estimated (optimized) vector.
+ * @param[in] (*get_f_error) pointer to the error function.
+ * @param[in] (*get_jacobian) pointer to the Jacobian matrix.
+ *
+ * @return required iteration number.
+ *
+ */
+uint8_t modified_gauss_newton(uint8_t f_length, uint8_t n,
+ vector_t x0_vec[n],
+ vector_t data_vec[f_length],
+ matrix_t eps, matrix_t fmin, uint8_t max_iter_num,
+ vector_t est_x_vec[n],
+ void (*get_f_error)(vector_t x0_vec[],
+ vector_t data_vec[],
+ vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x0_vec[],
+ matrix_t J[][n])
+ );
+
+#endif /* MODIFIED_GAUSS_NEWTON_H_ */
diff --git a/non_linear_algebra/optimization/levenberg_marquardt.c b/non_linear_algebra/optimization/levenberg_marquardt.c
new file mode 100644
index 0000000000000000000000000000000000000000..d1ce59df0521de8d2c176c0617ad51f90ffddb13
--- /dev/null
+++ b/non_linear_algebra/optimization/levenberg_marquardt.c
@@ -0,0 +1,239 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Implement the Levenberg--Marquardt (LVM) algorithm.
+ * @note This function is generally implemented.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * @author Naouar Guerchali
+ *
+ * @}
+ */
+
+#include <stdio.h>
+#include <math.h>
+#include <float.h>
+#include <inttypes.h>
+
+#include "levenberg_marquardt.h"
+#include "matrix.h"
+#include "vector.h"
+#include "solve.h"
+#include "utils.h"
+
+/**
+ * @brief Implements the correction-function of the Levenberg--Marquardt (LVM) algorithm.
+ * @details The user should provide pointers to the error and Jacobian functions.
+ *
+ * @note This function is generally implemented.
+ *
+ * @param[in] f_length length of the error functions vector.
+ * @param[in] n length of the start vector.
+ * @param[in] x_vec[] start vector.
+ * @param[in] data_vec[] data vector.
+ * @param[in] mu regularization parameter \f$ \mu \f$.
+ * @param[out] s[] correction vector.
+ * @param[in] (*get_f_error) pointer to the error function that calculates the matrix
+ * \f$ J_f^{T} J_{f} \f$.
+ * @param[in] (*get_jacobian) pointer to the Jacobian matrix.
+ *
+ *
+ * @return the parameter \f$ \rho_{\mu} \f$
+ */
+matrix_t opt_levenberg_marquardt_correction(uint8_t f_length, uint8_t n,
+ matrix_t x_vec[n],
+ matrix_t data_vec[f_length],
+ matrix_t mu,
+ matrix_t s[n],
+ void (*get_f_error)(
+ vector_t x_vec[],
+ vector_t data_vec[],
+ vector_t f_vec[]),
+ void (*get_jacobian)(
+ vector_t x_vec[],
+ matrix_t J[][n])
+ )
+{
+
+ matrix_t Fx[f_length];
+ matrix_t J[f_length][n];
+ vector_t f_vec[f_length];
+ matrix_t JTJ_mu2_I[n][n];
+ vector_t JTF[n];
+ matrix_t ro_mu = 0;
+ matrix_t Fx_square = 0;
+ matrix_t Fx_plus_s_square = 0;
+ matrix_t Fx_plus_J_mul_s_square = 0;
+ vector_t Fx_plus_J_mul_s[f_length];
+ vector_t F_x_plus_s[f_length];
+ vector_t x_plus_s[n];
+ matrix_t denom;
+
+ /*
+ * compute: -(J'J + mu^2*I)
+ */
+ //JT_J = J'*J
+ get_jacobian(x_vec, J);
+
+ //store J'J in JTJ_mu2_I
+ matrix_trans_mul_itself(f_length, n, J, JTJ_mu2_I);
+
+ //compute: J'J + mu^2*I
+ matrix_add_to_diag(n, JTJ_mu2_I, n, pow(mu, 2));
+
+ // -(J'*J+mu^2*eye(n)), is an (nxn) matrix
+ matrix_mul_scalar(n, n, JTJ_mu2_I, -1, JTJ_mu2_I);
+
+ //JT_f = J'*f
+ get_f_error(x_vec, data_vec, f_vec);
+ matrix_trans_mul_vec(f_length, n, J, f_length, f_vec, JTF);
+
+ //solve the equation: s = -(J'*J+mu^2*I)\(J'*Fx);
+ solve_householder(n, n, JTJ_mu2_I, JTF, s);
+
+ //f(x)
+ get_f_error(x_vec, data_vec, Fx);
+
+ //f(x)^2
+ Fx_square = vector_get_scalar_product(f_length, Fx, Fx);
+
+ //(x+s)
+ vector_add(n, x_vec, s, x_plus_s);
+
+ //f(x+s)
+ get_f_error(x_plus_s, data_vec, F_x_plus_s);
+
+ //f(x+s)^2
+ Fx_plus_s_square = vector_get_scalar_product(f_length, F_x_plus_s,
+ F_x_plus_s);
+
+ //compute: J(x)*s
+ matrix_mul_vec(f_length, n, J, s, Fx_plus_J_mul_s);
+
+ //compute: f(x) + J(x)*s
+ vector_add(f_length, Fx, Fx_plus_J_mul_s, Fx_plus_J_mul_s);
+
+ //compute: ||f(x) + J(x)*s||^2
+ Fx_plus_J_mul_s_square = vector_get_scalar_product(f_length,
+ Fx_plus_J_mul_s,
+ Fx_plus_J_mul_s);
+
+ denom = Fx_square - Fx_plus_J_mul_s_square;
+ if (denom != 0) {
+ ro_mu = (Fx_square - Fx_plus_s_square) / denom;
+ }
+ else {
+ puts("ro_mu is infinite !!!");
+ ro_mu = FLT_MAX;
+ }
+
+ return ro_mu;
+}
+
+uint8_t opt_levenberg_marquardt(uint8_t f_length, uint8_t n,
+ vector_t x0_vec[n],
+ vector_t data_vec[f_length],
+ matrix_t eps, matrix_t tau, matrix_t beta0,
+ matrix_t beta1,
+ uint8_t max_iter_num,
+ vector_t est_x_vec[n],
+ void (*get_f_error)(vector_t x0_vec[],
+ vector_t data_vec[],
+ vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x0_vec[],
+ matrix_t J[][n])
+ )
+{
+ matrix_t J[f_length][n];
+ vector_t JT_f[n];
+ matrix_t JTJ_mu2_I[n][n];
+ vector_t s[n];
+ vector_t f_vec[f_length];
+ matrix_t mu;
+ matrix_t ro_mu;
+ uint8_t it;
+
+ /*
+ * compute mu0 and -(J'J + mu0^2*I)
+ */
+ //JT_J = J'*J
+ get_jacobian(x0_vec, J);
+
+ //store J'J in JTJ_mu2_I
+ matrix_trans_mul_itself(f_length, n, J, JTJ_mu2_I);
+
+ // compute mu_0: mu = mu_0
+ mu = opt_levenberg_marquardt_get_mu0(n, tau, JTJ_mu2_I);
+
+ //compute: J'J + mu0^2*I
+ matrix_add_to_diag(n, JTJ_mu2_I, n, pow(mu, 2));
+
+ //-(J'*J+mu0^2*eye(n)) is a (n,n) matrix
+ matrix_mul_scalar(n, n, JTJ_mu2_I, -1, JTJ_mu2_I);
+
+ //compute JTF: JT_f = J'*f
+ get_f_error(x0_vec, data_vec, f_vec);
+ matrix_trans_mul_vec(f_length, n, J, f_length, f_vec, JT_f);
+
+ //solve the equation: s=-(J'*J+mu^2*eye(n))\(J'*F);
+ solve_householder(n, n, JTJ_mu2_I, JT_f, s);
+
+ vector_copy(n, x0_vec, est_x_vec);
+
+ it = 0;
+ while ((vector_get_norm2(n, s)
+ > eps * (1 + vector_get_norm2(n, est_x_vec)))
+ && (it < max_iter_num)) { //norm(s,2)>eps*(1+norm(x,2))
+
+ ro_mu = opt_levenberg_marquardt_correction(f_length, n,
+ est_x_vec, data_vec,
+ mu, s, get_f_error,
+ get_jacobian);
+
+ while (1) {
+ if (ro_mu <= beta0) {
+ mu = 2.0 * mu;
+ ro_mu = opt_levenberg_marquardt_correction(
+ f_length, n,
+ est_x_vec,
+ data_vec, mu, s, get_f_error,
+ get_jacobian);
+ }
+ else if (ro_mu >= beta1) {
+ mu = mu / 2.0;
+ break;
+ }
+ else {
+ break;
+ }
+ }
+ vector_add(n, est_x_vec, s, est_x_vec);
+ it = it + 1;
+ } //while
+
+ return it;
+}
+
+matrix_t opt_levenberg_marquardt_get_mu0(uint8_t n, matrix_t tau,
+ matrix_t JTJ[][n])
+{
+ matrix_t max_diag_JTJ = JTJ[0][0];
+
+ for (uint8_t i = 1; i < n; i++) {
+ if (JTJ[i][i] > max_diag_JTJ) {
+ max_diag_JTJ = JTJ[i][i];
+ }
+ }
+
+ return tau * max_diag_JTJ;
+}
diff --git a/non_linear_algebra/optimization/modified_gauss_newton.c b/non_linear_algebra/optimization/modified_gauss_newton.c
new file mode 100644
index 0000000000000000000000000000000000000000..37a40c10901fe81cfa7bdd6f640e82f514793f4c
--- /dev/null
+++ b/non_linear_algebra/optimization/modified_gauss_newton.c
@@ -0,0 +1,113 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Implement the Gauss--Newton algorithm.
+ * @note This function is generally implemented.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * @author Abdelmoumen Norrdine <a.norrdine@googlemail.com>
+ *
+ * @}
+ */
+
+#include <stdio.h>
+#include <math.h>
+
+#include "utils.h"
+#include "matrix.h"
+#include "vector.h"
+#include "moore_penrose_pseudo_inverse.h"
+
+//n is the size of x0 that is equal to the column number of J
+uint8_t modified_gauss_newton(uint8_t f_length, uint8_t n,
+ vector_t x0_vec[n],
+ vector_t data_vec[f_length],
+ matrix_t eps, matrix_t fmin, uint8_t max_iter_num,
+ vector_t est_x_vec[n],
+ void (*get_f_error)(vector_t x0_vec[],
+ vector_t data_vec[],
+ vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x0_vec[],
+ matrix_t J[][n])
+ )
+
+{
+ matrix_t J[f_length][n];
+ matrix_t JT_J[n][n];
+ matrix_t JT_f[n];
+ vector_t f_vec[f_length];
+ matrix_t f_error;
+ matrix_t pinv_JTJ_mat[n][n];
+ vector_t correction_vec[n];
+ vector_t x_vec[n];
+ vector_t next_x_vec[n];
+ matrix_t max_error, min_error;
+ matrix_t step;
+ uint8_t iter_num;
+
+ get_f_error(x0_vec, data_vec, f_vec);
+ f_error = vector_get_norm2(f_length, f_vec);
+ max_error = f_error;
+ min_error = max_error;
+ step = eps;
+
+ vector_copy(n, x0_vec, x_vec);
+ vector_copy(n, x0_vec, est_x_vec);
+ iter_num = 0;
+
+ while ((step >= eps) && (iter_num < max_iter_num) && (f_error > fmin)) {
+ /*
+ * Compute then correction terms & next x_vec values
+ */
+ //JT_J = J'*J
+ get_jacobian(x_vec, J);
+ matrix_trans_mul_itself(f_length, n, J, JT_J);
+
+ //JT_f = J'*f
+ get_f_error(x_vec, data_vec, f_vec);
+ matrix_trans_mul_vec(f_length, n, J, f_length, f_vec, JT_f);
+
+ //solve: J'J*s = -J'*f
+ moore_penrose_get_pinv(n, n, JT_J, pinv_JTJ_mat);
+ //s = (J'J)\J'*f
+ matrix_mul_vec(n, n, pinv_JTJ_mat, JT_f, correction_vec);
+
+ // x = x - s
+ vector_sub(n, x_vec, correction_vec, next_x_vec);
+
+ //next step
+ step = vector_get_euclidean_distance(n, x_vec, next_x_vec);
+
+ // x_vec = next_x_vec
+ vector_copy(n, next_x_vec, x_vec);
+
+ //error vector
+ get_f_error(x_vec, data_vec, f_vec);
+
+ f_error = vector_get_norm2(f_length, f_vec);
+
+ if (min_error > f_error) { //store the x_vec value with the minimum error in est_x_vec
+ vector_copy(n, x_vec, est_x_vec);
+ min_error = f_error; // update min_error
+ }
+
+ max_error = utils_max(f_error, max_error); // update max_error
+ iter_num++;
+ if ((max_error - min_error) > 10) {
+ break;
+ }
+
+ } //while
+
+ return iter_num;
+}
diff --git a/non_linear_algebra/solve_non_linear_equations/damped_newton_raphson.c b/non_linear_algebra/solve_non_linear_equations/damped_newton_raphson.c
new file mode 100644
index 0000000000000000000000000000000000000000..6833629a510932da9ad7fe7b7075620c74071363
--- /dev/null
+++ b/non_linear_algebra/solve_non_linear_equations/damped_newton_raphson.c
@@ -0,0 +1,138 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Implement the damped Newton--Raphson algorithm.
+ * @details The damped Newton--Raphson algorithm enables to solve
+ * multi-variant nonlinear equation systems.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <stdio.h>
+
+#include "moore_penrose_pseudo_inverse.h"
+#include "damped_newton_raphson.h"
+#include "vector.h"
+#include "matrix.h"
+
+//Multivariate Newton�s Method
+//n is the size of x0 that is equal to the column number of J
+uint8_t damped_newton_raphson(uint8_t f_length, uint8_t n, vector_t x0_arr[],
+ double min_lamda, double eps, uint8_t max_it_num,
+ vector_t est_x_arr[],
+ void (*get_non_lin_sys)(vector_t x_arr[],
+ vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x_arr[],
+ matrix_t J[][n]))
+{
+ uint8_t iter_num;
+ double lamda;
+ double lamda_c;
+ double damped_norm_x;
+ double damped_norm_prev_x;
+ vector_t prev_x_arr[n];
+ vector_t delta_x[n];
+ vector_t lamda_mul_delta_x[n];
+
+ damped_norm_x = get_damped_norm(f_length, n, x0_arr, get_non_lin_sys,
+ get_jacobian);
+ iter_num = 0;
+ vector_copy(n, x0_arr, prev_x_arr);
+ while ((damped_norm_x >= eps) && (iter_num < max_it_num)) {
+
+ get_delta_x(f_length, n, prev_x_arr, get_non_lin_sys,
+ get_jacobian, delta_x);
+ lamda = 1.0;
+
+ // x = x_k +lamda*s_k (lamda=1)
+ vector_add(n, prev_x_arr, delta_x, est_x_arr);
+
+ lamda_c = 1 - lamda / 4;
+ //|||f(x)|||_k
+ damped_norm_x = get_damped_norm(f_length, n, est_x_arr,
+ get_non_lin_sys, get_jacobian);
+
+ //|||f(prev_x)|||_k
+ damped_norm_prev_x = vector_get_norm2(n, delta_x);
+
+ while (damped_norm_x > (lamda_c * damped_norm_prev_x)) {
+ lamda = lamda / 2;
+ if (lamda >= min_lamda) {
+ //x = x_k +lamda*s_k
+ vector_scalar_mul(n, delta_x, lamda,
+ lamda_mul_delta_x);
+ vector_add(n, prev_x_arr, lamda_mul_delta_x,
+ est_x_arr);
+ }
+ else {
+ break;
+ }
+ // |||f(x)|||_k
+ damped_norm_x = get_damped_norm(f_length, n, est_x_arr,
+ get_non_lin_sys, get_jacobian);
+
+ } //while
+
+ if (lamda <= min_lamda) {
+ break;
+ }
+
+ //update prev_x
+ vector_copy(n, est_x_arr, prev_x_arr);
+
+ // |||f(x)|||_k
+ damped_norm_x = get_damped_norm(f_length, n, est_x_arr,
+ get_non_lin_sys, get_jacobian);
+ iter_num++;
+ }
+
+ return iter_num;
+}
+
+double get_damped_norm(uint8_t m, uint8_t n, vector_t x_arr[],
+ void (*get_non_lin_sys)(vector_t x_arr[],
+ vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n])
+ )
+{
+ double damp_norm = 0.0;
+ vector_t delta_x_arr[n];
+
+ get_delta_x(m, n, x_arr, get_non_lin_sys, get_jacobian, delta_x_arr);
+ damp_norm = vector_get_norm2(n, delta_x_arr);
+
+ return damp_norm;
+}
+
+//Computes the correction vector.
+void get_delta_x(uint8_t m, uint8_t n, vector_t x_arr[],
+ void (*get_non_lin_sys)(vector_t x_arr[], vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n]),
+ vector_t delta_x_arr[])
+{
+ matrix_t J[m][n];
+ matrix_t pinv_J[n][m];
+ vector_t f_vec[m];
+
+ // compute the Jacobian matrix for the vector x_arr.
+ get_jacobian(x_arr, J);
+ // compute the inverse of J (J^-1)
+ moore_penrose_get_pinv(m, n, J, pinv_J);
+ // compute the value of f(x)
+ get_non_lin_sys(x_arr, f_vec);
+ //delta_x = J1^-1*f(x)
+ matrix_mul_vec(n, m, pinv_J, f_vec, delta_x_arr);
+ vector_in_place_scalar_mul(n, delta_x_arr, -1);
+}
diff --git a/non_linear_algebra/solve_non_linear_equations/doc.txt b/non_linear_algebra/solve_non_linear_equations/doc.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a2d1bec0745f4226061e8e85eb858847bfddc43c
--- /dev/null
+++ b/non_linear_algebra/solve_non_linear_equations/doc.txt
@@ -0,0 +1,20 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser
+ * General Public License v2.1. See the file LICENSE in the top level
+ * directory for more details.
+ */
+
+/**
+ * @defgroup solve_non_linear_equations SOLVE_NON_LINEAR_EQUATIONS
+ * @ingroup non_linear_algebra
+ *
+ * @brief Enables to solve multi-variant nonlinear equation systems.
+ *
+ * @details The multi-variant nonlinear equation systems are solved using
+ * damped or the Newton--Raphson algorithms.
+ *
+ * @author Zakaria Kasmi
+ */
diff --git a/non_linear_algebra/solve_non_linear_equations/fsolve.c b/non_linear_algebra/solve_non_linear_equations/fsolve.c
new file mode 100644
index 0000000000000000000000000000000000000000..3b55379ec546270ac17b2aa8426d89086d739afe
--- /dev/null
+++ b/non_linear_algebra/solve_non_linear_equations/fsolve.c
@@ -0,0 +1,62 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Solve multi-variant nonlinear equation systems.
+ * @details The multi-variant nonlinear equation systems are solved using
+ * damped or the Newton--Raphson algorithms.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include "vector.h"
+#include "matrix.h"
+#include "fsolve.h"
+#include "newton_raphson.h"
+#include "damped_newton_raphson.h"
+
+uint8_t fsolve(uint8_t f_length, uint8_t x0_length, vector_t x0_arr[],
+ enum NON_LIN_ALGORITHM algo, vector_t est_x_arr[],
+ void (*get_non_lin_sys)(vector_t x_arr[], vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x_arr[], matrix_t J[][x0_length]))
+{
+ double tol = 1e-9;
+ uint8_t max_it_num = 77;
+ uint8_t iter_num = 0;
+
+ switch (algo) {
+ case Newton_Raphson:
+ iter_num = newton_raphson(f_length, x0_length, x0_arr, tol,
+ max_it_num,
+ est_x_arr, get_non_lin_sys, get_jacobian);
+ break;
+
+ case Damped_Newton_Raphson:
+ {
+ double min_lamda = 4.8828125e-04;
+ iter_num = damped_newton_raphson(f_length, x0_length, x0_arr,
+ min_lamda,
+ tol, max_it_num, est_x_arr,
+ get_non_lin_sys, get_jacobian);
+ break;
+ }
+
+ default:
+ iter_num = newton_raphson(f_length, x0_length, x0_arr, tol,
+ max_it_num,
+ est_x_arr, get_non_lin_sys, get_jacobian);
+ }
+
+ return iter_num;
+}
diff --git a/non_linear_algebra/solve_non_linear_equations/include/damped_newton_raphson.h b/non_linear_algebra/solve_non_linear_equations/include/damped_newton_raphson.h
new file mode 100644
index 0000000000000000000000000000000000000000..1fd9806e2f0b4a3e1a63e72903b5ae4ff1c4f45d
--- /dev/null
+++ b/non_linear_algebra/solve_non_linear_equations/include/damped_newton_raphson.h
@@ -0,0 +1,94 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup solve_non_linear_equations
+ * @{
+ *
+ * @file
+ * @brief Implement the damped Newton--Raphson algorithm.
+ * @details The damped Newton--Raphson algorithm enables to solve
+ * multi-variant nonlinear equation systems.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef DAMPED_NEWTON_RAPHSON_H_
+#define DAMPED_NEWTON_RAPHSON_H_
+
+#include "matrix.h"
+#include "vector.h"
+
+/**
+ * @brief Implements the damped Newton--Raphson algorithm.
+ * @details The user should provide pointers to non-linear equation systems and Jacobian functions.
+ * @note This function is generally implemented.
+ *
+ * @param[in] f_length length of the error functions vector.
+ * @param[in] n length of the start vector.
+ * @param[in] x0_arr[] start vector.
+ * @param[in] min_lamda minimal damping factor.
+ * @param[in] eps accuracy bound.
+ * @param[in] max_it_num maximal iteration number of the damped Newton--Raphson
+ * algorithm.
+ * @param[out] est_x_arr[] estimated (solution) vector.
+ * @param[in] (*get_non_lin_sys) pointer to non-linear equation systems.
+ * @param[in] (*get_jacobian) pointer to the Jacobian matrix.
+ *
+ * @return required iteration number.
+ *
+ */
+uint8_t damped_newton_raphson(uint8_t f_length, uint8_t n, vector_t x0_arr[],
+ double min_lamda, double eps, uint8_t max_it_num,
+ vector_t est_x_arr[],
+ void (*get_non_lin_sys)(vector_t x_arr[],
+ vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x_arr[],
+ matrix_t J[][n]));
+/**
+ * @brief Compute the norm of the damped Newton--Raphson algorithm.
+ * @details The user should provide pointers to non-linear equation systems and Jacobian functions.
+ * @note This function is generally implemented.
+ *
+ * @param[in] m number of the non-linear equations.
+ * @param[in] n length of the guess vector.
+ * @param[in] x_arr[] guess vector.
+ * @param[in] (*get_non_lin_sys) pointer to non-linear equation systems.
+ * @param[in] (*get_jacobian) pointer to the Jacobian matrix.
+ *
+ * @return norm of the damped Newton--Raphson algorithm.
+ *
+ */
+double get_damped_norm(uint8_t m, uint8_t n, vector_t x_arr[],
+ void (*get_non_lin_sys)(vector_t x_arr[],
+ vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n])
+ );
+
+/**
+ * @brief Compute the correction vector the damped Newton--Raphson algorithm.
+ * @details The user should provide pointers to non-linear equation systems and Jacobian functions.
+ * @note This function is generally implemented.
+ *
+ * @param[in] m number of the non-linear equations.
+ * @param[in] n length of the guess vector.
+ * @param[in] x_arr[] guess vector.
+ * @param[in] (*get_non_lin_sys) pointer to non-linear equation systems.
+ * @param[in] (*get_jacobian) pointer to the Jacobian matrix.
+ * @param[in, out] delta_x_arr[] the correction vector (term).
+ *
+ */
+void get_delta_x(uint8_t m, uint8_t n, vector_t x_arr[],
+ void (*get_non_lin_sys)(vector_t x_arr[], vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n]),
+ vector_t delta_x_arr[]);
+
+#endif /* DAMPED_NEWTON_RAPHSON_H_ */
diff --git a/non_linear_algebra/solve_non_linear_equations/include/fsolve.h b/non_linear_algebra/solve_non_linear_equations/include/fsolve.h
new file mode 100644
index 0000000000000000000000000000000000000000..35f27cb8b1b84f38aff9514f65c87e79a6159c1e
--- /dev/null
+++ b/non_linear_algebra/solve_non_linear_equations/include/fsolve.h
@@ -0,0 +1,58 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup solve_non_linear_equations
+ * @{
+ *
+ * @file
+ * @brief Solve multi-variant nonlinear equation systems.
+ * @details The multi-variant nonlinear equation systems are solved using
+ * damped or the Newton--Raphson algorithms.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef FSOLVE_H_
+#define FSOLVE_H_
+
+/**
+ * Possible algorithms to solve multi-variant nonlinear equation systems.
+ */
+enum NON_LIN_ALGORITHM {
+ Newton_Raphson, /**< Newton--Raphson algorithm */
+ Damped_Newton_Raphson /**< Damped Newton--Raphson algorithm */
+};
+
+/**
+ * @brief Solve systems of multi-variant nonlinear equations.
+ * @details The user should provide pointers to non-linear equation systems and Jacobian functions.
+ * The user can choose between the damped or the Newton--Raphson algorithms.
+ *
+ * @note This function is generally implemented.
+ *
+ * @param[in] f_length length of the error functions vector.
+ * @param[in] x0_length length of the start vector.
+ * @param[in] x0_arr[] start vector.
+ * @param[in] algo damped or the Newton--Raphson algorithm.
+ * @param[in, out] est_x_arr[] estimated (solution) vector.
+ * @param[in] (*get_non_lin_sys) pointer to non-linear equation systems.
+ * @param[in] (*get_jacobian) pointer to the Jacobian matrix.
+ *
+ * @return required iteration number.
+ *
+ */
+uint8_t fsolve(uint8_t f_length, uint8_t x0_length, vector_t x0_arr[],
+ enum NON_LIN_ALGORITHM algo, vector_t est_x_arr[],
+ void (*get_non_lin_sys)(vector_t x_arr[], vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x_arr[], matrix_t J[][x0_length]));
+
+#endif /* FSOLVE_H_ */
diff --git a/non_linear_algebra/solve_non_linear_equations/include/newton_raphson.h b/non_linear_algebra/solve_non_linear_equations/include/newton_raphson.h
new file mode 100644
index 0000000000000000000000000000000000000000..2a5edf98fb7cec0e133cc11b329b6c7ad511d18d
--- /dev/null
+++ b/non_linear_algebra/solve_non_linear_equations/include/newton_raphson.h
@@ -0,0 +1,53 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup solve_non_linear_equations
+ * @{
+ *
+ * @file
+ * @brief Implement the Newton--Raphson algorithm.
+ * @details The Newton--Raphson algorithm enables to solve
+ * multi-variant nonlinear equation systems.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+#ifndef NEWTON_RAPHSON_H_
+#define NEWTON_RAPHSON_H_
+
+#include <inttypes.h>
+
+#include "vector.h"
+#include "matrix.h"
+
+/**
+ * @brief Implements the Newton--Raphson algorithm.
+ * @details The user should provide pointers to non-linear equation systems and Jacobian functions.
+ * @note This function is generally implemented.
+ *
+ * @param[in] f_length length of the error functions vector.
+ * @param[in] n length of the start vector.
+ * @param[in] x0_arr[] start vector.
+ * @param[in] eps accuracy bound.
+ * @param[in] max_it_num maximal iteration number of the Newton--Raphson algorithm.
+ * @param[out] est_x_arr[] estimated (solution) vector.
+ * @param[in] (*get_non_lin_sys) pointer to non-linear equation systems.
+ * @param[in] (*get_jacobian) pointer to the Jacobian matrix.
+ *
+ * @return required iteration number.
+ *
+ */
+uint8_t newton_raphson(uint8_t f_length, uint8_t n, vector_t x0_arr[],
+ double eps, uint8_t max_it_num, vector_t est_x_arr[],
+ void (*get_non_lin_sys)(vector_t x_arr[], vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n]));
+
+#endif /* NEWTON_RAPHSON_H_ */
diff --git a/non_linear_algebra/solve_non_linear_equations/newton_raphson.c b/non_linear_algebra/solve_non_linear_equations/newton_raphson.c
new file mode 100644
index 0000000000000000000000000000000000000000..fd07c893755792ec9f7b6ca4c7bd6b698839c123
--- /dev/null
+++ b/non_linear_algebra/solve_non_linear_equations/newton_raphson.c
@@ -0,0 +1,65 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Implement the Newton--Raphson algorithm.
+ * @details The Newton--Raphson algorithm enables to solve
+ * multi-variant nonlinear equation systems.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include "vector.h"
+#include "matrix.h"
+#include "moore_penrose_pseudo_inverse.h"
+
+//n is the size of x0 that is equal to the column number of J
+uint8_t newton_raphson(uint8_t f_length, uint8_t n, vector_t x0_arr[],
+ double eps, uint8_t max_it_num, vector_t est_x_arr[],
+ void (*get_non_lin_sys)(vector_t x_arr[],
+ vector_t f_vec[]),
+ void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n]))
+{
+ uint8_t iter_num;
+ double step;
+ vector_t prev_x_arr[n];
+ matrix_t J[f_length][n];
+ matrix_t pinv_J[n][f_length];
+ vector_t delta_x[n];
+ vector_t f_vec[f_length];
+
+ vector_copy(n, x0_arr, prev_x_arr);
+ step = eps;
+ iter_num = 0;
+
+ while ((step >= eps) && (iter_num < max_it_num)) {
+ // compute the Jacobian matrix for the prev_x_arr.
+ get_jacobian(prev_x_arr, J);
+ // compute the inverse of J (J^-1)
+ moore_penrose_get_pinv(f_length, n, J, pinv_J);
+ // compute the value of f(x)
+ get_non_lin_sys(prev_x_arr, f_vec);
+ //delta_x = J1^-1*f(x)
+ matrix_mul_vec(n, f_length, pinv_J, f_vec, delta_x);
+ //improve x: x = x - delta_x
+ vector_sub(n, prev_x_arr, delta_x, est_x_arr);
+ //next step
+ step = vector_get_euclidean_distance(n, est_x_arr, prev_x_arr);
+ //update prev_x
+ vector_copy(n, est_x_arr, prev_x_arr);
+ iter_num++;
+ }
+
+ return iter_num;
+}