diff --git a/Cholesky/CMakeLists.txt b/Cholesky/CMakeLists.txt
new file mode 100644
index 0000000..e7fcacc
--- /dev/null
+++ b/Cholesky/CMakeLists.txt
@@ -0,0 +1,10 @@
+include_directories(${CMAKE_CURRENT_SOURCE_DIR})
+
+# The next lines represent a list of all the executables in the folder (chapter)
+# There is typically one executable per exercise.
+# As you solve exercises, you need to add them here.
+# The syntax is add_executable(executable_name source_file_name)
+
+
+add_executable(main main.cpp Cholesky.cpp Functions.cpp)
+               
\ No newline at end of file
diff --git a/Cholesky/Cholesky.cpp b/Cholesky/Cholesky.cpp
new file mode 100644
index 0000000..9828037
--- /dev/null
+++ b/Cholesky/Cholesky.cpp
@@ -0,0 +1,39 @@
+#include "Cholesky.hpp"
+
+//Default constructor//
+Cholesky::Cholesky(double **m_A, double *m_b, int m_size) {
+    A=m_A;
+    b=m_b;
+    size=m_size;
+    Cholesky_matrices Structure=Cholesky_L(A,b,size);
+    L=Structure.L_matrix;
+    L_T=Structure.L_T_matrix;
+    b_star=Structure.b_star_vector;
+}
+
+//Copy constructor//
+Cholesky::Cholesky(const Cholesky &otherCholesky) {
+    A=otherCholesky.A;
+    b=otherCholesky.b;
+    b_star=otherCholesky.b_star;
+    size=otherCholesky.size;
+    L=otherCholesky.L;
+    L_T=otherCholesky.L_T;
+}
+
+double** Cholesky::get_L() {
+    return L;
+}
+
+double** Cholesky::get_L_T() {
+    return L_T;
+}
+
+double* Cholesky::get_b_star(){
+    return b_star;
+}
+
+double* Cholesky::Solve() {
+    double* z=Up(L_T,b_star,size);
+    return z;
+}
\ No newline at end of file
diff --git a/Cholesky/Cholesky.hpp b/Cholesky/Cholesky.hpp
new file mode 100644
index 0000000..1a9f4f2
--- /dev/null
+++ b/Cholesky/Cholesky.hpp
@@ -0,0 +1,24 @@
+#ifndef CHOLESKYHEADERREF
+#define CHOLESKYHEADERREF
+
+#include "Functions.hpp"
+
+class Cholesky{
+private:
+    int size;
+    double** A;
+    double* b;
+    double* b_star;
+    double** L;
+    double** L_T;
+
+public:
+    Cholesky(const Cholesky& otherCholesky);
+    Cholesky(double** m_A,double* m_b, int m_size);
+    double** get_L();
+    double** get_L_T();
+    double* get_b_star();
+    double* Solve();
+};
+
+#endif
\ No newline at end of file
diff --git a/Cholesky/Functions.cpp b/Cholesky/Functions.cpp
new file mode 100644
index 0000000..ea76ffe
--- /dev/null
+++ b/Cholesky/Functions.cpp
@@ -0,0 +1,292 @@
+#include <iostream>
+#include <cassert>
+#include <cmath>
+#include "Functions.hpp"
+#include "Cholesky.hpp"
+
+double** AllocateMemory(int rows, int cols){
+    /*Allocates memory for matrix A*/
+    double** A;
+    A=new double* [rows];
+    for (int i=0;i<rows;i++){
+        A[i]=new double[cols];
+    }
+    return A;
+}
+
+double** Initialize(int size){
+    double** A=AllocateMemory(size,size);
+    for (int i=0;i<size;i++) {
+        for (int j = 0; j < size; j++) {
+            A[i][j] = 0;
+        }
+    }
+    for (int i=0;i<size;i++){
+        A[i][i]=1;
+    }
+    return A;
+}
+
+double** Copy(double** A,int size){
+    double** B=AllocateMemory(size,size);
+    for (int i=0;i<size;i++){
+        for (int j=0;j<size;j++){
+            B[i][j]=A[i][j];
+        }
+    }
+    return B;
+}
+
+double* Copy(double* b,int size){
+    double* b_star=AllocateMemory(size);
+    for (int i=0;i<size;i++){
+            b_star[i]=b[i];
+        }
+    return b_star;
+}
+
+void SubstractRows (double** A,double scalar,int row_1, int row_2, int size){
+    for (int j=0;j<size;j++){
+        A[row_1][j]-=A[row_2][j]*scalar;
+    }
+}
+
+void SubstractRows (double* u,double scalar,int row_1, int row_2){
+        u[row_1]-=u[row_2]*scalar;
+}
+
+double Normalize (double** A,int row, int size){
+    double scalar=A[row][row];
+    for (int j=0;j<size;j++){
+        A[row][j]/=scalar;
+    }
+    return scalar;
+}
+
+LU_matrices Pivot(double** A,double* b,int size) {
+        double** L=Initialize(size);
+        double** U=Copy(A,size);
+        double* b_star=Copy(b,size);
+
+        LU_matrices output;
+        for (int i=0;i<size-1;i++){
+            for (int j=i+1;j<size;j++){
+                L[j][i]=U[j][i]/U[i][i];
+                SubstractRows(U,U[j][i]/U[i][i],j,i,size);
+                SubstractRows(b_star,L[j][i],j,i);
+            }
+        }
+        output.U_matrix=U;
+        output.L_matrix=L;
+        output.b_star_vector=b_star;
+        return output;
+}
+
+
+double* AllocateMemory(int rows){
+    /*Allocates memory for vector u*/
+    double* u=new double [rows];
+    return u;
+}
+
+void FreeMemory(double** A,int rows){
+    /*Frees memory of matrix A*/
+    for (int i=0;i<rows;i++){
+        delete[] A[i];
+    }
+    delete[] A;
+}
+
+void FreeMemory(double* u){
+    /*Frees memory of vector u*/
+    delete[] u;
+}
+
+int FindMax (double** A,int col,int size){
+    int index=0;
+    for (int i=0;i<size;i++){
+        if (A[i][col]>A[index][col]){
+            index=i;
+        }
+    }
+    return index;
+}
+
+void SwapRows (double** A,int last_index, int new_index, int size){
+    double* temp=AllocateMemory(size);
+    for (int i=0;i<size;i++){
+        temp[i]=A[last_index][i];
+        A[last_index][i]=A[new_index][i];
+        A[new_index][i]=temp[i];
+    }
+    FreeMemory(temp);
+}
+
+
+
+void TerminalPrint(double** A,int rows, int cols){
+    /*Prints to terminal values of a matrix A*/
+    for (int i=0;i<rows;i++){
+        for (int j=0;j<cols;j++){
+            std::cout<<A[i][j]<<" ";
+        }
+        std::cout<<std::endl;
+    }
+    std::cout<<std::endl;
+}
+
+void TerminalPrint(double* u,int size){
+    /*Prints to terminal values of a vector u*/
+    for (int i=0;i<size;i++) {
+        std::cout << u[i]<<std::endl;
+    }
+    std::cout<<std::endl;
+    }
+
+double** Multiply(double** A,double** B,int rowsA, int colsA, int rowsB, int colsB){
+    /*Multiplies matrix A[rowsA][colsA] by matrix B[rowsB][colsB]*/
+    assert(rowsB==colsA);
+    double** C=AllocateMemory(rowsA,colsB);
+    for (int i=0;i<rowsA;i++){
+        for (int j=0;j<colsB;j++){
+            for (int k=0;k<rowsB;k++){
+                C[i][j]+=A[i][k]*B[k][j];
+            }
+        }
+    }
+    return C;
+}
+
+double* Multiply(double** A,double* u,int rowsA, int colsA, int rowsU){
+    /*Multiplies matrix A[rowsA][colsA] by vector U[rowsU]*/
+    assert(rowsU==colsA);
+    double* v=AllocateMemory(rowsA);
+    for (int i=0;i<rowsA;i++){
+        for (int j=0;j<colsA;j++){
+            v[i]+=A[i][j]*u[j];
+        }
+    }
+    return v;
+}
+
+double* Multiply(double* u,double** A,int rowsA, int colsA, int rowsU){
+    /*Multiplies vector U[rowsU] by matrix A[rowsA][colsA]*/
+    assert(rowsU==rowsA);
+    double* v=AllocateMemory(colsA);
+    for (int i=0;i<colsA;i++){
+        for (int j=0;j<rowsA;j++){
+            v[i]+=A[j][i]*u[j];
+        }
+    }
+    return v;
+}
+
+double** Multiply(double** A,double scalar, int rowsA, int colsA){
+    /*Multiplies matrix A[rowsA][colsA] by scalar*/
+    double** C=AllocateMemory(rowsA,colsA);
+    for (int i=0;i<rowsA;i++){
+        for (int j=0;j<colsA;j++){
+            C[i][j]=A[i][j]*scalar;
+        }
+    }
+    return C;
+}
+
+double* Multiply(double* u,double scalar, int rowsU){
+    /*Multiplies vector u[rowsU] by scalar*/
+    double* v=AllocateMemory(rowsU);
+    for (int i=0;i<rowsU;i++){
+      v[i]=u[i]*scalar;
+    }
+    return v;
+}
+
+void ReduceMatrix(double** A,int size,int col, double** &A_hat){
+    /*Resizing of matrix A for determinant computation*/
+    for (int row = 1; row < size; ++row) {
+        for (int k = 0; k < col; ++k) {
+            A_hat[row - 1][k] = A[row][k];
+        }
+        for (int k = col + 1; k < size; ++k) {
+            A_hat[row - 1][k - 1] = A[row][k];
+        }
+    }
+}
+
+double Determinant(double** A,int size){
+    /*Computation of determinant : direct formula for 2x2, recursive otherwise*/
+
+    if (size==2){
+        return A[0][0]*A[1][1]-A[1][0]*A[0][1];
+
+    } else {
+        double det=0;
+        for (int j = 0; j < size; ++j) {
+            double** A_hat=AllocateMemory(size-1,size-1);
+            ReduceMatrix(A,size,j,A_hat);
+            det +=A[0][j]*std::pow(-1.0,j)*Determinant(A_hat,size-1);
+           FreeMemory(A_hat,size-1);
+        }
+        return det;
+    }
+}
+
+double* Up(double** U, double* b_star, int size){
+    double* x=AllocateMemory(size);
+    for (int i=2;i>-1;i--){
+        x[i]+=b_star[i];
+        for (int j=i+1;j<size;j++){
+            x[i]-=U[i][j]*x[j];
+        }
+        x[i]/=U[i][i];
+    }
+return x;
+}
+
+double* Down(double** L, double* b_star, int size){
+    double* x=AllocateMemory(size);
+    for (int i=0;i<size;i++){
+        x[i]+=b_star[i];
+        for (int j=i-1;j>=0;j--){
+            x[i]-=L[i][j]*x[j];
+        }
+        x[i]/=L[i][i];
+    }
+    return x;
+}
+
+double** Transpose(double** A,int size){
+    double** temp=Copy(A,size);
+    for (int i=0;i<size;i++){
+        for (int j=0;j<size;j++){
+            temp[i][j]=A[j][i];
+        }
+    }
+    return temp;
+}
+
+Cholesky_matrices Cholesky_L(double** A,double* b,int size){
+    double** L=AllocateMemory(size,size);
+    double* b_star=Copy(b,size);
+    Cholesky_matrices output;
+
+    for (int i=0;i<size;i++){
+        for (int k=0;k<i+1;k++){
+            double sum=0;
+                for (int j=0;j<k;j++){
+                    sum+=(L[i][j]*L[k][j]);
+                }
+                if (i==k){
+                    L[i][k]=sqrt(A[i][i]-sum);
+                } else {
+                    L[i][k]=(1.0/L[k][k])*(A[i][k]-sum);
+                }
+            }
+        }
+    output.L_matrix=L;
+    output.L_T_matrix=Transpose(L,size);
+    output.b_star_vector=Down(L,b,size);
+    return output;
+
+}
+
diff --git a/Cholesky/Functions.hpp b/Cholesky/Functions.hpp
new file mode 100644
index 0000000..9c2a514
--- /dev/null
+++ b/Cholesky/Functions.hpp
@@ -0,0 +1,43 @@
+#ifndef FUNCTIONSHEADERREF
+#define FUNCTIONSHEADERREF
+
+
+class LU_matrices{
+public:
+    double** U_matrix;
+    double** L_matrix;
+    double* b_star_vector;
+};
+
+class Cholesky_matrices{
+public:
+    double** L_matrix;
+    double** L_T_matrix;
+    double* b_star_vector;
+};
+
+double** Initialize(int size);
+double** Copy(double** A,int size);
+double** AllocateMemory(int rows, int cols);
+void SubstractRows (double** A,int scalar,int row_1, int row_2, int size);
+double Normalize (double** A,int row, int size);
+LU_matrices Pivot(double** A,double* b,int size);
+double* AllocateMemory(int rows);
+void FreeMemory(double** A,int rows);
+void FreeMemory(double* u);
+int FindMax (double** A,int col,int size);
+void SwapRows (double** A,int last_index, int new_index, int size);
+void TerminalPrint(double** A,int rows, int cols);
+void TerminalPrint(double* u,int size);
+double** Multiply(double** A,double** B,int rowsA, int colsA, int rowsB, int colsB);
+double* Multiply(double** A,double* u,int rowsA, int colsA, int rowsU);
+double* Multiply(double* u,double** A,int rowsA, int colsA, int rowsU);
+double** Multiply(double** A,double scalar, int rowsA, int colsA);
+double* Multiply(double* u,double scalar, int rowsU);
+void ReduceMatrix(double** A,int size,int col, double** &A_hat);
+double Determinant(double** A,int size);
+double** Generate_L(double** L,LU_matrices out,int size);
+double* Up(double** U, double* b_star, int size);
+Cholesky_matrices Cholesky_L(double** A,double* b,int size);
+
+#endif
diff --git a/Cholesky/main.cpp b/Cholesky/main.cpp
new file mode 100644
index 0000000..17884cd
--- /dev/null
+++ b/Cholesky/main.cpp
@@ -0,0 +1,40 @@
+#include <iostream>
+#include <cassert>
+#include <cmath>
+#include "Functions.hpp"
+#include "Cholesky.hpp"
+
+
+int main(int argc, char *argv[]){
+    int size=3;
+    double** A=AllocateMemory(size,size);
+    double* b=AllocateMemory(size);
+
+    A[0][0]=25;A[0][1]=15;A[0][2]=-5;
+    A[1][0]=15;A[1][1]=18;A[1][2]=0;
+    A[2][0]=-5;A[2][1]=0;A[2][2]=11;
+
+    b[0]=1;b[1]=6;b[2]=4;
+    std::cout<<"Matrix A:"<<std::endl;
+    TerminalPrint(A,size,size);
+    Cholesky linsys(A,b,size);
+    std::cout<<"Matrix L:"<<std::endl;
+    TerminalPrint(linsys.get_L(),size,size);
+    std::cout<<"Matrix L_T:"<<std::endl;
+    TerminalPrint(linsys.get_L_T(),size,size);
+    double** A_back=Multiply(linsys.get_L(),linsys.get_L_T(),size,size,size,size);
+    std::cout<<"Matrix L * Matrix L_T:"<<std::endl;
+    TerminalPrint(A_back,size,size);
+    std::cout<<"Vector b:"<<std::endl;
+    TerminalPrint(b,size);
+    std::cout<<"Vector z:"<<std::endl;
+    TerminalPrint(linsys.get_b_star(),size);
+    std::cout<<"Solution of system Ax=b:"<<std::endl;
+    TerminalPrint(linsys.Solve(),size);
+
+    return 0;
+
+
+
+
+}