diff --git a/doc/notebooks/Gaussian_Quadrature.ipynb b/doc/notebooks/Gaussian_Quadrature.ipynb
new file mode 100644
index 000000000..6c82e8620
--- /dev/null
+++ b/doc/notebooks/Gaussian_Quadrature.ipynb
@@ -0,0 +1,170 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Gaussian Quadrature\n",
+ "\n",
+ "The general idea of Gaussian quadrature is to find, in an interval, a set of points and weights such that the integral of a function in the interval will be approximate to the sum of the weighted function values at the points. That is:\n",
+ "\n",
+ "\\begin{equation*}\n",
+ " \\int_{a}^{b}f(x) \\simeq \\sum_{i=1}^{n}w_{i}f(x_{i})\n",
+ " \\tag{1}\n",
+ "\\end{equation*} \n",
+ " \n",
+ "\n",
+ "This method is designed to give highly accurate results for $n$ points given that the function to integrate can be approximated accurately as a multiple of a polynomial of grade $m \\leq 2n-1$ by a weight function $w(x)$.\n",
+ "\n",
+ "\n",
+ "\\begin{equation}\n",
+ " f(x)\\simeq L^{m}w(x)\n",
+ " \\tag{2}\n",
+ "\\end{equation}\n",
+ "\n",
+ "The way to select the points and weights can be understood as follows:\n",
+ "\n",
+ "Given a set of $n$ polynomials $\\{p^{i}(x)\\}$ orthogonal in the interval $[a,b]$ with respect to a weight function $w(x)$:\n",
+ "\n",
+ "\\begin{equation}\n",
+ " \\int_{a}^{b} p^{i}(x) p^{j}(x) w(x) = C_{i} \\delta_{ij}\n",
+ " \\tag{3}\n",
+ "\\end{equation}\n",
+ "\n",
+ "Where $\\delta_{ij}$ is the Kronecker delta:\n",
+ "\n",
+ "\\begin{equation}\n",
+ " \\delta_{ij} = \\begin{cases}\n",
+ " 0 & \\text{if } i \\neq j, \\\\\n",
+ " 1 & \\text{if } i = j.\n",
+ " \\end{cases}\n",
+ " \\tag{4}\n",
+ "\\end{equation}\n",
+ "\n",
+ "Any polynomial of order $m\\leq2n-1$ $(L^{m}(x))$ divided by $p^{n}(x)$ will give a quotient $Q(x)$ and remainder $R(x)$ of order at most $n-1$. Then $L^{m}(x)$ can be expressed as:\n",
+ "\\begin{equation}\n",
+ " L^{m}(x) = Q(x)p^{n}(x) + R(x) \n",
+ " \\tag{5}\n",
+ "\\end{equation}\n",
+ "\n",
+ "Also, as the remainder is of order lower or equal to $n-1$, it can be expressed as a linear combination of the orthogonal polynomials $\\{p^{j}(x)\\}$ for $j=0\\ldots n-1$\n",
+ "\n",
+ "\\begin{equation}\n",
+ " R(x) = \\sum_{j=0}^{n-1}c_{j}p^{j}(x)\n",
+ " \\tag{6}\n",
+ "\\end{equation}\n",
+ "\n",
+ "Substituting [5](#eqn:poly_div) in [2](#eqn:poly_approx) gives:\n",
+ "\n",
+ "\\begin{equation}\n",
+ " f(x)\\simeq [Q(x)p^{n}(x) + R(x)]w(x)\n",
+ " \\tag{7}\n",
+ "\\end{equation}\n",
+ "\n",
+ "Substituting [6](#eqn:remaind_expanx):\n",
+ "\n",
+ "\\begin{equation}\n",
+ " f(x)\\simeq [Q(x)p^{n}(x) + \\sum_{j=0}^{n-1}c_{j}p^{j}(x)]w(x)\n",
+ " \\tag{8}\n",
+ "\\end{equation}\n",
+ "\n",
+ "\n",
+ "Then the approximate function integral can be expressed as:\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\\begin{align}\n",
+ " \\int_{a}^{b}f(x) &\\simeq \\int_{a}^{b}[Q(x)p^{n}(x) + \\sum_{j=0}^{n-1}c_{j}p^{j}(x)]w(x)dx \\tag{9} \\\\\n",
+ " &\\simeq \\int_{a}^{b}Q(x)p^{n}(x)w(x) dx + \\int_{a}^{b}\\sum_{j=0}^{n-1}c_{j}p^{j}(x)w(x)dx \\tag{10} \\\\ \n",
+ " &\\simeq \\int_{a}^{b}Q(x)p^{n}(x)w(x) dx + \\sum_{j=0}^{n-1}c_{j}\\int_{a}^{b}p^{j}(x)w(x)dx \\tag{11} \t\t\n",
+ "\\end{align} \n",
+ "\n",
+ "where the first integral on the right side is 0 because $p^{n}(x)$ is orthogonal to any polynomial of lower order in the $[a,b]$ interval with respect to the weight function $w(x)$:\n",
+ "\n",
+ "\\begin{equation}\n",
+ " \\int_{a}^{b}Q(x)p^{n}(x)w(x) dx = 0\n",
+ " \\tag{12}\n",
+ "\\end{equation}\n",
+ "\n",
+ "Also substituting [8](#eqn:fcn_aprox_poly_div_exp) in [1](#eqn:num_int):\n",
+ "\n",
+ "\\begin{align}\n",
+ " \\int_{a}^{b}f(x) \\simeq \\sum_{i=1}^{n}w_{i}f(x_{i}) &= \\sum_{i=1}^{n}w_{i} \\left[Q(x_{i})p^{n}(x_{i}) + \\sum_{j=0}^{n-1}c_{j}p^{j}(x_{i})\\right]w(x_{i}) \\tag{13}\\\\\n",
+ " &= \\sum_{i=1}^{n}w_{i} Q(x_{i})p^{n}(x_{i}) + \\sum_{i=1}^{n}w_{i}\\sum_{j=0}^{n-1}c_{j}p^{j}(x_{i})w(x_{i}) \\tag{14}\\\\\n",
+ " &= \\sum_{i=1}^{n}w_{i} Q(x_{i})p^{n}(x_{i}) + \\sum_{j=0}^{n-1}c_{j}\\sum_{i=1}^{n}w_{i}p^{j}(x_{i})w(x_{i}) \\tag{15}\n",
+ "\\end{align} \n",
+ "\n",
+ "To make numeric the integral of the function as close to the exact value as possible the $\\{x_{i},w_{i}\\}$ values are selected such that the sums in [15](#eq:num_int_dec) are equal to the integrals in [11](#eqn:aprox_int_exp_reord). The $\\{x_{i}\\}$ vales are then selected to be the roots of $p^{n}(x)$. In this way the leftmost sum is zero (the exact value of the integral [12](#eqn:orth_quot)) for any possible set of weights $\\{w_{i}\\}$.\n",
+ "\\begin{equation}\n",
+ " \\sum_{i=1}^{n}w_{i} Q(x_{i})p^{n}(x_{i}) = 0 \\tag{16}\n",
+ "\\end{equation}\n",
+ "\n",
+ "The corresponding weights are then selected such that for all $n-1$ polynomials $p^{j}(x)$, the sum give their exact integral value:\n",
+ "\\begin{equation}\n",
+ " \\sum_{i=1}^{n}w_{i}p^{j}(x_{i})w(x_{i}) = \\int_{a}^{b}p^{j}(x)w(x)dx \\tag{17}\n",
+ "\\end{equation}\n",
+ "\n",
+ "This condition guarantees that the rightmost sum in [15](#eq:num_int_dec) gives the exact value of the rightmost integral in [11](#eqn:aprox_int_exp_reord). \n",
+ "\n",
+ "The weights that fulfill [17](#eqn:int_sum_weights) can be then found by solving the following linear system.\n",
+ "\n",
+ "\\begin{equation}\n",
+ " \\begin{bmatrix}\n",
+ " p^{0}(x_{1})w(x_{1})&p^{0}(x_{2})w(x_{2})&p^{0}(x_{3})&\\dots &p^{0}(x_{n})w(x_{n})\\\\\n",
+ " p^{1}(x_{1})w(x_{1})&p^{1}(x_{2})w(x_{2})&p^{1}(x_{3})&\\dots &p^{1}(x_{n})w(x_{n})\\\\\n",
+ " p^{2}(x_{1})w(x_{1})&p^{2}(x_{2})w(x_{2})&p^{2}(x_{3})&\\dots &p^{2}(x_{n})w(x_{n})\\\\\n",
+ " \\vdots &\\vdots &\\vdots &\\ddots &\\vdots \\\\\n",
+ " p^{(n-1)}(x_{1})w(x_{1})&p^{(n-1)}(x_{2})w(x_{2})&p^{(n-1)}(x_{3})&\\dots &p^{(n-1)}(x_{n})w(x_{n})\\\\\n",
+ " \\end{bmatrix}\n",
+ " \\cdot \n",
+ " \\begin{bmatrix}\n",
+ " w_{0}\\\\w_{1}\\\\w_{2}\\\\\\vdots \\\\w_{n}\n",
+ " \\end{bmatrix}\n",
+ " =\n",
+ " \\begin{bmatrix}\n",
+ " \\int_{a}^{b} p^{0}(x)w(x)\\\\\n",
+ " \\int_{a}^{b} p^{1}(x)w(x)\\\\\n",
+ " \\int_{a}^{b} p^{2}(x)w(x)\\\\\n",
+ " \\vdots \\\\\n",
+ " \\int_{a}^{b} p^{(n-1)}(x)w(x)\\\\\n",
+ " \\end{bmatrix}\n",
+ " \\tag{18}\n",
+ "\\end{equation}\n",
+ "\n",
+ "Given a family of orthogonal polynomials with respect to a weight function in a determined interval, this procedure allows obtaining a set of $n$ optimal points $\\{x_{i}\\}$ and their corresponding weights $\\{w_{i}\\}$ for integrating a function in this interval. The function needs to be approximateable as a polynomial of grade at most $2n-1$. Also, if the function can be represented exactly as a polynomial of grade at most $2n-1$, the numerical integration is exact.\n",
+ "\n",
+ "There are several Gauss quadrature rules for different choices of orthogonal polynomial sets, intervals and weight functions. Grid offers several one-dimensional quadratures for one-dimensional integration such as:\n",
+ "\n",
+ "\t\n",
+ "- Gauss-Chebyshev\n",
+ "\n",
+ "- Gauss-Legendre\n",
+ "\n",
+ "- Gauss-Laguerre"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "grid_env",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.12"
+ },
+ "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}