diff --git a/MLPP/.DS_Store b/MLPP/.DS_Store
new file mode 100644
index 0000000..342bf4f
Binary files /dev/null and b/MLPP/.DS_Store differ
diff --git a/MLPP/NumericalAnalysis/NumericalAnalysis.cpp b/MLPP/NumericalAnalysis/NumericalAnalysis.cpp
index c4b812c..52f6332 100644
--- a/MLPP/NumericalAnalysis/NumericalAnalysis.cpp
+++ b/MLPP/NumericalAnalysis/NumericalAnalysis.cpp
@@ -14,6 +14,12 @@ namespace MLPP{
         return (function(x + eps) - function(x)) / eps; // This is just the formal def. of the derivative.
     }
 
+    
+    double NumericalAnalysis::numDiff_2(double(*function)(double), double x){ 
+        double eps = 1e-5;
+        return (function(x + eps) -  2 * function(x) + function(x - eps)) / (eps * eps);
+    }
+
     double NumericalAnalysis::numDiff(double(*function)(std::vector<double>), std::vector<double> x, int axis){
         // For multivariable function analysis. 
         // This will be used for calculating Jacobian vectors. 
@@ -25,9 +31,26 @@ namespace MLPP{
         return (function(x_eps) - function(x)) / eps; 
     }
 
-    double NumericalAnalysis::newtonRaphsonMethod(double(*function)(double), double x_0, double epoch){
+    double NumericalAnalysis::numDiff_2(double(*function)(std::vector<double>), std::vector<double> x, int axis1, int axis2){
+        //For Hessians. 
+        double eps = 1e-5;
+
+        std::vector<double> x_pp = x;
+        x_pp[axis1] += eps; 
+        x_pp[axis2] += eps; 
+
+        std::vector<double> x_np = x;
+        x_np[axis2] += eps; 
+            
+        std::vector<double> x_pn = x;
+        x_pn[axis1] += eps;
+
+        return (function(x_pp) - function(x_np) - function(x_pn) + function(x))/(eps * eps);
+    }
+
+    double NumericalAnalysis::newtonRaphsonMethod(double(*function)(double), double x_0, double epoch_num){
         double x = x_0;
-        for(int i = 0; i < epoch; i++){
+        for(int i = 0; i < epoch_num; i++){
             x = x - function(x)/numDiff(function, x);
         }
         return x;
@@ -41,4 +64,17 @@ namespace MLPP{
         }
         return jacobian;
     }
+    std::vector<std::vector<double>> NumericalAnalysis::hessian(double(*function)(std::vector<double>), std::vector<double> x){
+        std::vector<std::vector<double>> hessian; 
+        hessian.resize(x.size());
+        for(int i = 0; i < hessian.size(); i++){
+            hessian[i].resize(x.size());
+        }
+        for(int i = 0; i < hessian.size(); i++){
+            for(int j = 0; j < hessian[i].size(); j++){
+                hessian[i][j] = numDiff_2(function, x, i, j);
+            }
+        }
+        return hessian;
+    }
 }
\ No newline at end of file
diff --git a/MLPP/NumericalAnalysis/NumericalAnalysis.hpp b/MLPP/NumericalAnalysis/NumericalAnalysis.hpp
index f6e6977..ee3313f 100644
--- a/MLPP/NumericalAnalysis/NumericalAnalysis.hpp
+++ b/MLPP/NumericalAnalysis/NumericalAnalysis.hpp
@@ -16,10 +16,15 @@ namespace MLPP{
             the future.
             */
             double numDiff(double(*function)(double), double x);
-            double numDiff(double(*function)(std::vector<double>), std::vector<double> x, int axis);
-            double newtonRaphsonMethod(double(*function)(double), double x_0, double epoch);
+            double numDiff_2(double(*function)(double), double x);
 
-            std::vector<double> jacobian(double(*function)(std::vector<double>), std::vector<double>);
+            double numDiff(double(*function)(std::vector<double>), std::vector<double> x, int axis);
+            double numDiff_2(double(*function)(std::vector<double>), std::vector<double> x, int axis1, int axis2);
+
+            double newtonRaphsonMethod(double(*function)(double), double x_0, double epoch_num);
+
+            std::vector<double> jacobian(double(*function)(std::vector<double>), std::vector<double> x); // Indeed, for functions with scalar outputs the Jacobians will be vectors.
+            std::vector<std::vector<double>> hessian(double(*function)(std::vector<double>), std::vector<double> x);
         
     };
 }
diff --git a/a.out b/a.out
new file mode 100755
index 0000000..02c5923
Binary files /dev/null and b/a.out differ
diff --git a/main.cpp b/main.cpp
index a5e2424..a847ee0 100644
--- a/main.cpp
+++ b/main.cpp
@@ -50,14 +50,43 @@ using namespace MLPP;
 
 
 double f(double x){
-    return x*x*x + 2*x - 2;
+    return x*x*x + 2*x - 2; 
 }
+/*
+    y = x^3 + 2x - 2
+    y' = 3x^2 + 2
+    y'' = 6x
+    y''(2) = 12
+*/
+
+// double f_mv(std::vector<double> x){
+//     return x[0] * x[0] + x[0] * x[1] * x[1] + x[1] + 5; 
+// }
+
+/* 
+    Where x, y = x[0], x[1], this function is defined as:
+    f(x, y) = x^2 + xy^2 + y + 5
+    ∂f/∂x = 2x + 2y
+    ∂^2f/∂x∂y = 2
+*/
 
 double f_mv(std::vector<double> x){
-    return x[0] * x[0] + x[1] * x[1] + x[1] + 5; 
-    // Where x,y=x[0],x[1], this function is defined as:
-    // f(x,y) = x^2 + y^2 + y + 5
+    return x[0] * x[0] * x[0] + x[0] + x[1] * x[1] * x[1] * x[0] + x[2] * x[2] * x[1];
 }
+/*
+    Where x, y = x[0], x[1], this function is defined as:
+    f(x, y) = x^4 + xy^3 + yz^2
+
+    ∂f/∂x = 3x^2 + 3y^2 
+    ∂^2f/∂x^2 = 6x 
+
+    ∂f/∂z = 2zy
+    ∂^2f/∂z^2 = 2z
+
+    ∂f/∂y = 3xy^2
+    ∂^2f/∂y∂x = 3y^2
+*/
+
 
 int main() {
 
@@ -475,12 +504,18 @@ int main() {
     // Checks for numerical analysis class.
     NumericalAnalysis numAn;
 
-    std::cout << numAn.numDiff(&f, 1) << std::endl;
-    std::cout << numAn.newtonRaphsonMethod(&f, 1, 1000) << std::endl;
+    // std::cout << numAn.numDiff(&f, 1) << std::endl;
+    // std::cout << numAn.newtonRaphsonMethod(&f, 1, 1000) << std::endl;
 
-    std::cout << numAn.numDiff(&f_mv, {1, 1}, 1) << std::endl; // Derivative w.r.t. x.
+    // std::cout << numAn.numDiff(&f_mv, {1, 1}, 1) << std::endl; // Derivative w.r.t. x.
 
-    alg.printVector(numAn.jacobian(&f_mv, {1, 1}));
+    // alg.printVector(numAn.jacobian(&f_mv, {1, 1}));
+
+    //std::cout << numAn.numDiff_2(&f, 2) << std::endl;
+
+    std::cout << numAn.numDiff_2(&f_mv, {2, 2, 500}, 2, 2) << std::endl;
+    std::cout << "Our Hessian." << std::endl;
+    alg.printMatrix(numAn.hessian(&f_mv, {2, 2, 500}));
 
     return 0;
 }