diff --git a/Gauss-Elim/GEtest b/Gauss-Elim/GEtest
new file mode 100644
index 0000000..40d2f79
Binary files /dev/null and b/Gauss-Elim/GEtest differ
diff --git a/Gauss-Elim/GEtest.cpp b/Gauss-Elim/GEtest.cpp
new file mode 100644
index 0000000..79eb07f
--- /dev/null
+++ b/Gauss-Elim/GEtest.cpp
@@ -0,0 +1,41 @@
+// GEtest.cpp
+// Author: Isaac Stone
+// Language: C++
+//
+// Synopsis: This program simply tests the return
+// value of the gaussElim() function described in
+// "gaussElim.h". The mat_gen() function generates
+// an augmented matrix that is stored into A. Both 
+// the matrix and the solution to Ax=b are printed
+// to the screen
+
+#include <stdlib.h>
+#include <stdio.h>
+#include "matrix.h"
+#include "vectorCode.h"
+#include "mat_gen.h"
+#include "gaussElim.h"
+
+using namespace std;
+
+int main(  int argc, char **argv ) {
+
+  Matrix A;
+  vector<double> a_sols;
+  
+  A = mat_gen();
+  a_sols.resize( A.size(), 0 );
+  
+  // Perform Gussian Elimination with 
+  // back substitution
+  a_sols = gaussElim( A );
+
+  // Print augmented matrix A
+  printMatrix(A);
+  printf( "\n" );
+
+  // Printf solutions to Ax=b
+  printVector( a_sols );
+
+  return 0;
+}
diff --git a/Gauss-Elim/README.md b/Gauss-Elim/README.md
new file mode 100644
index 0000000..1233e34
--- /dev/null
+++ b/Gauss-Elim/README.md
@@ -0,0 +1,26 @@
+### Gaussian Elimination and Back-Substution Algorithm - gaussElim.h
+
+#### Author: Isaac Stone, University of Tennessee Undergrad of Computer Science
+
+### gaussElim.h
+gaussElim.h - Included is gaussElim.h, which implements the Gaussian Elimination and back-substituion method for an augmented matrix. The functionality is fairly simple. The function gaussElim() can be called by any c++ program that includes "gaussElim.h". It takes as its only argument a Matrix ( described in matrix.cpp ) and returns a Vector ( described in vectorCode.cpp ) that is the solution to Ax=b. If no solution exists, or only the trivial solution exists, the program prints an error and exits. 
+
+#### Other files worth noting:
+
+mat_gen.h - Generates an augmented matrix using pseudorandom values. It returns a Matrix and can be included via "mat_gen.h". For more, see commments in file.
+
+GEtest.cpp - Used to test gaussElim.h with various inputs from mat_gen.h.
+
+makefile - If working on a linux system, make sure all files are in working directory and type "make" to compile.
+ 
+iterativeSolvers.h, iterativeSolver.cpp, matrix.h, matrix.cpp, vectorCode.h, VectorCode.cpp -
+All these files are dependencies for Matrix and Vector data structures needed for gaussElim.h and mat_gen.h to function as-is. They are written and developed by Christopher Johnson and obtained via GitHub at:
+
+### https://github.com/Christopher42/computational-mathematics
+
+where extensive documentation can be found, along with examples and testing information.
+
+Project was compiled on a linux machine using g++ -std=11 and was produced under the GNU GeneralPublic License, version 3. Feel free to modify, include, and improve gaussElim.h and mat_gen.h under the GGPL. Thank you. 
+
+
+Acknowledgements to Christopher Johnson, the author of the computational-mathematics repository on GitHhub.
diff --git a/Gauss-Elim/gaussElim.cpp b/Gauss-Elim/gaussElim.cpp
new file mode 100644
index 0000000..9102572
--- /dev/null
+++ b/Gauss-Elim/gaussElim.cpp
@@ -0,0 +1,89 @@
+// GEtest.cpp
+// Author: Isaac Stone
+// Language: C++
+//
+// Synopsis: This .cpp file implements Guassin Elimination with
+// back-substitution for an augmented matrix. The function
+// gaussElim() takes a Matrix ( matrix.h ) as an argument and returns a
+// solution to Ax=b in the form of a Vector ( vectorCode.h ). If no solution or 
+// only the trivial solution exists, a message is printed 
+// to the screen and the program exits.
+
+#include <vector>
+#include <cstdio>
+#include <cstdlib>
+#include "matrix.h"
+#include "vectorCode.h"
+
+using namespace std;
+
+// Elimination step
+void Elim( double x, Vector &a, Vector &b ) {
+	
+  for ( unsigned i = 0; i < a.size(); i++ )
+	  a[i] = a[i] - ( x * b[i] );
+}
+
+Vector gaussElim( Matrix A ) {
+	
+  int i, j, p, n;
+  double sum, xn, x;
+  Vector tmp, sols;
+
+  x = 0.0;
+  p = 0;
+  n = A.size();
+
+  // Guassian Elimination
+  for ( i = 0; i < (n - 1) ; i++ )
+  {
+	p = i;
+	while( A[p][i] == 0 )
+	{
+	  p++;
+	  if ( p >  n )
+	  {
+		printf( "No solution or trivial solution\n" );
+		exit(1);
+	  }
+	}
+	
+	// Swap rows
+	if ( p != i ) {
+	  tmp = A[p];
+	  A[p] = A[i];
+	  A[i] = tmp;
+	}
+		
+	for ( j = i+1; j < n; j++ )
+	{
+	  x = A[j][i] / A[i][i];
+	  Elim( x, A[j], A[i] );
+	}
+  }
+
+  // Test for no solution
+  if ( A[ n - 1 ][ n - 1 ] == 0 )
+  {
+    printf( "No solution or trivial solution\n" );
+	exit(1);
+  }
+ 
+  // Perform back-substitution
+  sols.resize( n, 0 );
+  
+  xn = A[n-1][n] / A[n-1][n-1];
+  sols[n-1] = xn;
+  
+  for ( int i = n-1; i >= 0; i-- )
+  {
+	sum = 0;
+	for ( j = i+1; j < n+1; j++ ) {
+	  sum += ( A[i][j] * sols[j] );
+	}
+	xn = ( A[i][n] - sum ) / A[i][i];
+	sols[i] = xn;
+  }
+
+  return sols;
+}
diff --git a/Gauss-Elim/gaussElim.h b/Gauss-Elim/gaussElim.h
new file mode 100644
index 0000000..d7b720e
--- /dev/null
+++ b/Gauss-Elim/gaussElim.h
@@ -0,0 +1,19 @@
+// GEtest.cpp
+// Author: Isaac Stone
+// Language: C++
+//
+// Synopsis: This file includes headers gaussElim.cpp to implement
+// Guassin Elimination with back-substitution for an augmented matrix. The function
+// gaussElim() takes a Matrix ( matrix.h ) as an argument and returns a
+// solution to Ax=b in the form of a Vector ( vectorCode.h ). If no solution or 
+// only the trivial solution exists, a message is printed 
+// to the screen and the program exits.
+//
+// Check "gaassElim.cpp" for more information.
+
+#include <vector>
+#include <cstdlib>
+#include "matrix.h"
+#include "vectorCode.h"
+
+Vector gaussElim( Matrix A );
diff --git a/Gauss-Elim/iterativeSolvers.cpp b/Gauss-Elim/iterativeSolvers.cpp
new file mode 100644
index 0000000..df66758
--- /dev/null
+++ b/Gauss-Elim/iterativeSolvers.cpp
@@ -0,0 +1,280 @@
+#include "matrix.h"
+#include "vectorCode.h"
+#include "iterativeSolvers.h"
+#include <cmath>
+
+int jacobi(Matrix const &A, Vector const &b, Vector & x, int maxIter, double tol)
+{
+	int n = A.size();
+	if (maxIter == 0) //makes default maxIter=n
+		maxIter = n;
+	double err = 10*tol;
+	Vector x_old = x;
+	int k=0;
+
+	while (k++<maxIter and err>tol*tol)
+	{
+		// x = D^-1((L+U)x-b)
+		for (int i=0;i<n;++i)
+		{
+			double sum = 0.0;
+			for (int j=0;j<i;++j)
+				sum+=A[i][j]*x_old[j];
+			for (int j=i+1;j<n;++j)
+				sum+=A[i][j]*x_old[j];
+			x[i]=(b[i]-sum)/A[i][i];
+		}
+
+		// err = error^2 = L2err(x,x_old)^2
+		err = 0.0;
+		for (int i=0;i<n;++i)
+			err += (x[i]-x_old[i])*(x[i]-x_old[i]);
+		x_old = x;
+	}
+	return k;
+}
+
+int gaussSeidel (Matrix const &A, Vector const &b, Vector & x, int maxIter, double tol)
+{
+	int n = A.size();
+	if (maxIter == 0)
+		maxIter = n;
+	double err = 10*tol;
+	Vector x_old = x;
+	int k=0;
+
+	while (k++<maxIter and err>tol*tol)
+	{
+		for (int i=0;i<n;++i)
+		{
+			double sum = 0;
+			for (int j=0;j<i;++j)
+				sum+=A[i][j]*x[j];
+			for (int j=i+1;j<n;++j)
+				sum+=A[i][j]*x[j];
+			x[i]=(b[i]-sum)/A[i][i];
+		}
+
+		// err = error^2 = L2err(x,x_old)^2
+		err = 0.0;
+		for (int i=0;i<n;++i)
+			err += (x[i]-x_old[i])*(x[i]-x_old[i]);
+		x_old = x;
+	}
+	return k;
+}
+
+int steepestDescent (Matrix const &A, Vector const &b, Vector & x, int maxIter, double tol)
+{
+	int n = A.size();
+	int k=0;
+	if (maxIter == 0)
+		maxIter = n;
+	//r_0 = b - Ax
+	Vector r = matrixVectorProduct(A,x);
+	for (int i=0;i<n;++i)
+		r[i] = b[i]-r[i];
+	//delta = ||x||_2^2
+	double delta = dotProduct(r,r);
+	//b_delta = ||b||_2^2
+	double b_delta = dotProduct(b,b);
+	Vector p = r;
+
+	while (k++<maxIter and delta>tol*tol*b_delta)
+	{
+		//s_k = A*p_k
+		Vector s(n,0);
+		for (int i=0;i<n;++i)
+			for (int j=0;j<n;++j)
+				s[i]+=A[i][j]*p[j];
+		//alpha_k = delta_k / (r_k \dot s_k)
+		double alpha = 0;
+		for (int i=0;i<n;++i)
+			alpha += r[i]*s[i];
+		alpha = delta/alpha;
+		//x_k+1 = x_k + alpha_k*p_k
+		for (int i=0;i<n;++i)
+			x[i] += alpha*p[i];
+		//r_k+1 = r_k - alpha_k*s_k
+		for (int i=0;i<n;++i)
+			r[i] -= alpha*s[i];
+		//delta_k+1 = ||r_k+1||_2^2
+		delta = 0;
+		for (int i=0;i<n;++i)
+			delta+=r[i]*r[i];
+		p=r;
+	}
+	return k;
+}
+
+int conjugateGradient (Matrix const &A, Vector const &b, Vector & x, int maxIter, double tol)
+{
+	int n = A.size();
+	if (maxIter == 0)
+		maxIter = n;
+	//r_0 = b - Ax
+	Vector x1(n,0);
+	Vector r = matrixVectorProduct(A,x);
+	for (int i=0;i<n;++i)
+		r[i] = b[i]-r[i];
+	//delta = ||r||_2^2
+	double delta = dotProduct(r,r);
+	//b_delta = ||b||_2^2
+	double b_delta = dotProduct(b,b);
+	
+	int k=0;
+	Vector p = r;
+
+	while (k++<maxIter and delta>tol*tol*b_delta)
+	{
+		//s_k = A*p_k
+		Vector s(n,0);
+		for (int i=0;i<n;++i)
+			for (int j=0;j<n;++j)
+				s[i]+=A[i][j]*p[j];
+		//alpha_k = delta_k / (p_k \dot s_k)
+		double alpha = 0;
+		for (int i=0;i<n;++i)
+			alpha += p[i]*s[i];
+		alpha = delta/alpha;
+		//x_k+1 = x_k + alpha_k*p_k
+		//r_k+1 = r_k - alpha_k*s_k
+		for (int i=0;i<n;++i)
+		{
+			x1[i] = x[i] + alpha*p[i];
+			r[i] = r[i] - alpha*s[i];
+		}
+			
+		//delta_k+1 = ||r_k+1||_2^2
+		double delta_new = 0;
+		for (int i=0;i<n;++i)
+			delta+=r[i]*r[i];
+		//p_k+1 = r_k+1 + (delta_k+1/delta_k)p_k
+		for (int i=0;i<n;++i)
+			p[i]=r[i]+(delta_new/delta)*p[i];
+		delta = delta_new;
+		x = x1;
+	}
+	return k;
+}
+
+//alt code for conjGradient with all the function calls
+// int conjugateGradient(Matrix const &A, Vector const &b,
+//         Vector &x0, int maxIter, double tol)
+// {   
+//     int n = b.size();
+//     if (maxIter == 0)
+//     	maxIter = n;
+//     Vector x1(n,0);
+//     Vector r0(n,0);
+   
+//     Vector Ax0 = matrixVectorProduct(A, x0);
+   
+//     for(int i = 0; i < n; i ++){
+//         r0[i] = b[i] - Ax0[i];
+//     }
+   
+//     double delta0 = dotProduct(r0,r0);
+//     double b_delta = dotProduct(b,b);
+   
+//     int k = 0;
+    
+//     Vector p0 = r0;
+   
+//     while(delta0 > tol*tol*b_delta and k < maxIter){
+       
+//         Vector s = matrixVectorProduct(A, p0);
+       
+//         double alpha = delta0/dotProduct(p0,s);
+     
+//         for(int i = 0; i < n; i++){
+//             x1[i] = x0[i] + alpha*p0[i];
+//             r0[i] = r0[i] - alpha*s[i];
+//         }
+       
+//         double delta1 = dotProduct(r0, r0);
+       
+//         for(int i = 0; i < n; i++){
+//             p0[i] = r0[i] + (delta1/delta0)*p0[i];
+//         }
+        
+//         k++;
+        
+//         delta0 = delta1;
+//         x0 = x1;
+        
+//     }
+//     return k;
+// }
+
+/*
+double powerMethod (Matrix const & A, Vector &x, int maxIter, double tol)
+{
+	int n = A.size();
+	Vector y = matrixVectorProduct(A,x);
+	double lambda = 0.0;
+	double err = 10*tol;
+	int k = 0;
+	while (k++<maxIter and err > tol)
+	{
+		// x =  y * (1/||y||);
+		double y_norm = 0.0;
+		for (int i=0; i<n; i++)
+			y_norm+=y[i]*y[i];
+		y_norm = sqrt(y_norm);
+		for(int i=0;i<n;i++)
+			x[i] = y[i]*(1/y_norm);
+		
+		// Vector s = A*x;
+		Vector s(n,0);
+		for (int i=0;i<n;++i)
+			for (int j=0;j<n;++j)
+				s[i]+=A[i][j]*x[j];
+
+		// lambda_new = dotProduct(x,s);
+		double lambda_new = 0.0;
+		for(int i=0;i<n;++i)
+			lambda_new+=x[i]*s[i];
+
+		err = abs(lambda - lambda_new);
+		lambda = lambda_new;
+	}
+	return lambda;
+}
+
+double inversePowerMethod (Matrix const & A, Vector &x, int maxIter, double tol)
+{
+	int n = A.size();
+	//s = A_inv*x;
+	Vector y(n,.5);
+	jacobi(A,x,y,40);
+	double lambda = 0.0;
+	double err = 10*tol;
+	int k = 0;
+	while (k++<maxIter and err > tol)
+	{
+		// x =  y * / ||y||;
+		// double y_norm = vectorNormL2(y);
+		// x = vectorScale(y,1.0/y_norm);
+		double y_norm = 0.0;
+		for (int i=0; i<n; i++)
+			y_norm+=y[i]*y[i];
+		y_norm = sqrt(y_norm);
+		for(int i=0;i<n;i++)
+			x[i] = y[i]*(1/y_norm);
+		
+		// Vector s = A_inv*x;
+		Vector s(n,.5);
+		jacobi(A,x,s,40);
+
+		// lambda_new = dotProduct(x,s);
+		double lambda_new = 0.0;
+		for(int i=0;i<n;++i)
+			lambda_new+=x[i]*s[i];
+
+		err = abs(lambda - lambda_new);
+		lambda = lambda_new;
+	}
+	return (1/lambda);
+}
+*/
diff --git a/Gauss-Elim/iterativeSolvers.h b/Gauss-Elim/iterativeSolvers.h
new file mode 100644
index 0000000..eaece8d
--- /dev/null
+++ b/Gauss-Elim/iterativeSolvers.h
@@ -0,0 +1,17 @@
+#ifndef ITERATIVE_SOLVERS_H
+#define ITERATIVE_SOLVERS_H
+
+#include <cmath>
+#include <vector>
+
+typedef std::vector<double> Vector;
+typedef std::vector<std::vector<double>> Matrix;
+
+int jacobi(Matrix const &A, Vector const &b, Vector & x, int maxIter=0, double tol=pow(10.0,-8.0));
+int gaussSeidel (Matrix const &A, Vector const &b, Vector & x, int maxIter=0, double tol=pow(10.0,-8.0));
+int steepestDescent (Matrix const &A, Vector const &b, Vector & x, int maxIter=0, double tol=pow(10.0,-8.0));
+int conjugateGradient (Matrix const &A, Vector const &b, Vector & x, int maxIter=0, double tol=pow(10.0,-10.0));
+
+//double powerMethod (Matrix const & A, Vector &x, int maxIter=100, double tol=pow(10.0,-10.0));
+//double inversePowerMethod (Matrix const & A, Vector &x, int maxIter=100, double tol=pow(10.0,-10.0));
+#endif
diff --git a/Gauss-Elim/makefile b/Gauss-Elim/makefile
new file mode 100644
index 0000000..d63123c
--- /dev/null
+++ b/Gauss-Elim/makefile
@@ -0,0 +1,4 @@
+
+
+make: 
+	g++ -std=c++11 -g GEtest.cpp gaussElim.cpp iterativeSolvers.cpp vectorCode.cpp matrix.cpp -o GEtest 
diff --git a/Gauss-Elim/mat_gen.h b/Gauss-Elim/mat_gen.h
new file mode 100644
index 0000000..5dac2dd
--- /dev/null
+++ b/Gauss-Elim/mat_gen.h
@@ -0,0 +1,51 @@
+// GEtest.cpp
+// Author: Isaac Stone
+// Language: C++
+//
+// Synopsis: This .h implements mat_gen() which  generates a pseudorandom augmented
+// matrix ( it may be solvable or not ) and returns as a Matrix ( "matrix.h" ).
+// Size of the matrix may be changed by modifying i and j limiting values within
+// for loops.
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <vector>
+#include "matrix.h"
+#include <random>
+#include <chrono>
+
+using namespace std;
+
+// set up random number generator
+default_random_engine dre (chrono::steady_clock::now().time_since_epoch().count());
+
+int random ( int lim ) {
+  uniform_int_distribution<int> uid {0,lim};
+  return uid(dre);
+}
+
+Matrix mat_gen() {
+
+  int i, j;
+  double x;
+  Matrix M;
+  vector <double> row;
+
+  // seed generator
+  srand( time(NULL) );
+
+  M.clear();
+
+  for (  i = 0; i < 5; i++ ) {
+
+    row.clear();
+	for ( j = 0; j < 6; j++ ) {
+	  x = random(20);
+	  row.push_back( x );
+	}
+    M.push_back( row );
+  }
+
+  return M;
+}
diff --git a/Gauss-Elim/matrix.cpp b/Gauss-Elim/matrix.cpp
new file mode 100644
index 0000000..b95d105
--- /dev/null
+++ b/Gauss-Elim/matrix.cpp
@@ -0,0 +1,116 @@
+#include "matrix.h"
+#include <cmath>
+#include <iostream>
+
+double matrixNormL1 (Matrix m)
+{
+	double maxColumnSum = 0;
+	for (int j=0;j<m[0].size();j++)
+	{
+		double sum = 0;
+		for (int i=0;i<m.size();i++)
+			sum+= std::abs(m[i][j]);
+		if (sum > maxColumnSum)
+			maxColumnSum = sum;
+	}
+	return maxColumnSum;
+}
+
+
+double matrixNormLInf (Matrix m)
+{
+	double maxRowSum = 0;
+	for (int i=0;i<m.size();i++)
+	{
+		double sum = 0;
+		for (int j=0;j<m[0].size();j++)
+			sum+= std::abs(m[i][j]);
+		if (sum > maxRowSum)
+			maxRowSum = sum;
+	}
+	return maxRowSum;
+}
+
+void printMatrix (Matrix m)
+{
+	for (int i=0;i<m.size();i++)
+	{
+		std::cout << ((i==0) ? '[' : ' ');
+		std::cout << "[" << m[i][0];
+		for (int j=1; j<m[i].size(); j++)
+			std::cout << "," << m[i][j];
+		std::cout << "]" << ((i==m.size()-1) ? ']' : ',') << std::endl;
+	}
+}
+
+Matrix identityMatrix(int n)
+{
+	Matrix ident(n, Vector(n,0));
+	for (int i=0;i<n;++i)
+		ident[i][i] = 1;
+	return ident;
+}
+
+Matrix matrixAdd(Matrix a, Matrix b)
+{
+	int n = a.size();
+	int m = a[0].size();
+	Matrix sol (m, Vector(n,0));
+
+	for (int i=0;i<n;++i)
+		for (int j=0;j<m;++j)
+			sol[i][j] = a[i][j] + b[i][j];
+
+	return sol;
+}
+
+Matrix matrixSub(Matrix a, Matrix b)
+{
+	int n = a.size();
+	int m = a[0].size();
+	Matrix sol (m, Vector(n,0));
+
+	for (int i=0;i<n;++i)
+		for (int j=0;j<m;++j)
+			sol[i][j] = a[i][j] - b[i][j];
+
+	return sol;
+}
+
+Matrix matrixScale(Matrix a, double s)
+{
+	Matrix sol = a;
+
+	for (int i=0;i<a.size();++i)
+		for (int j=0;j<a[0].size();j++)
+			sol[i][j]*= s;
+
+	return sol;
+}
+
+Vector matrixVectorProduct(Matrix a, Vector v)
+{
+	int n = a.size();
+	int m = a[0].size();
+	Vector sol(m,0);
+
+	for (int i=0;i<n;++i)
+		for (int j=0;j<m;++j)
+			sol[i]+=a[i][j]*v[j];
+	return sol;
+}
+
+Matrix matrixMatrixProduct(Matrix a, Matrix b)
+{
+	int n = a.size();
+	int m = b.size();
+	int r = b[0].size();
+
+	Matrix sol(n, Vector(r,0));
+
+	for(int i=0;i<n;++i)
+		for(int j=0;j<r;++j)
+			for(int k=0;k<m;++k)
+				sol[i][j]+=a[i][k]*b[k][j];
+	return sol;
+}
\ No newline at end of file
diff --git a/Gauss-Elim/matrix.h b/Gauss-Elim/matrix.h
new file mode 100644
index 0000000..7c5d7aa
--- /dev/null
+++ b/Gauss-Elim/matrix.h
@@ -0,0 +1,21 @@
+#ifndef MATRIX_H
+#define MATRIX_H
+
+#include <vector>
+
+typedef std::vector<double> Vector;
+typedef std::vector<std::vector<double>> Matrix;
+
+double matrixNormL1 (Matrix m);
+double matrixNormLInf (Matrix m);
+
+void printMatrix (Matrix m);
+Matrix identityMatrix(int n);
+
+Matrix matrixAdd(Matrix a, Matrix b);
+Matrix matrixSub(Matrix a, Matrix b);
+Matrix matrixScale(Matrix a, double s);
+Vector matrixVectorProduct(Matrix a, Vector v);
+Matrix matrixMatrixProduct(Matrix a, Matrix b);
+
+#endif
diff --git a/Gauss-Elim/test.h b/Gauss-Elim/test.h
new file mode 100644
index 0000000..1b89ec1
--- /dev/null
+++ b/Gauss-Elim/test.h
@@ -0,0 +1,90 @@
+#include <vector>
+#include <cstdlib>
+#include <iostream>
+
+using namespace std;
+
+void Swap( vector <double> &a, vector <double> &b ) {
+	
+	vector <double> tmp;
+	tmp = a;
+	a = b;
+	b = tmp;
+}
+
+void Elim( double x, vector <double> &a, vector <double> &b ) {
+	
+		for ( unsigned i = 0; i < a.size(); i++ )
+			a[i] = a[i] - ( x * b[i] );
+}
+
+int gaussElim( vector< vector <double> > &m, vector <double> &sols ) {
+	
+	unsigned int i, j, p, n, iterations;
+	double sum, xn, x;
+
+	x = 0.0;
+	p = 0;
+	iterations = 0;
+	n = m.size();
+	/*
+	for ( unsigned i = 0; i < m.size(); i++ )
+	{
+		for ( unsigned j = 0; j < m[0].size(); j++ )
+		{
+			cout << m[i][j] << " ";
+		}
+		cout << endl;
+	}
+	cout << endl << endl; 
+	*/
+
+	for ( i = 0; i < n - 1 ; i++ )
+	{
+	  iterations++;
+
+	  p = i;
+	  while( m[p][i] == 0 )
+	  {
+		p++;
+		if ( p > n )
+		{
+		  cout << "No solution or trivial solution" << endl;
+		  return 0;
+		}
+	  }
+	  
+	  if ( p != i )
+		Swap( m[p], m[i] );
+		  
+	  for ( j = i+1; j < n; j++ )
+	  {
+		x = m[j][i] / m[i][i];
+		Elim( x, m[j], m[i] );
+	  }
+	}
+	
+	if ( m[ n - 1 ][ n - 1 ] == 0 )
+	{
+	  cout << "No solution or trivial solution" << endl;
+		return 0;
+	}
+	
+	xn = m[n-1][n] / m[n-1][n-1];
+	sols.resize( n, 0 );
+	sols[n-1] = xn;
+	
+	for ( i = n-1; i >= 0; i-- )
+	{
+	  sum = 0;
+	  for ( j = i+1; j < n+1; j++ )
+		sum += m[i][j] * sols[j];
+	  xn = ( m[i][n] - sum ) / m[i][i];
+	  sols[i] = xn;
+	}
+
+	for ( i = 0; i < sols.size(); i++ )
+	  cout << - sols[i] << " ";
+	 
+	return iterations;
+}
diff --git a/Gauss-Elim/vectorCode.cpp b/Gauss-Elim/vectorCode.cpp
new file mode 100644
index 0000000..1b8049c
--- /dev/null
+++ b/Gauss-Elim/vectorCode.cpp
@@ -0,0 +1,76 @@
+#include "vectorCode.h"
+#include <cmath>
+#include <iostream>
+
+double vectorNormL1 (Vector v)
+{
+	double sum = 0.0;
+	for (int i=0; i<v.size(); i++)
+		sum+=std::abs(v[i]);
+	return sum;
+}
+
+double vectorNormL2 (Vector v)
+{
+	double sum_of_squares = 0.0;
+	for (int i=0; i<v.size(); i++)
+		sum_of_squares+=v[i]*v[i];
+	return sqrt(sum_of_squares);
+}
+
+double vectorNormLInf (Vector v)
+{
+	double maxElement = v[0];
+	for(int i=1; i<v.size(); i++)
+		if (v[i]>maxElement)
+			maxElement=v[i];
+	return maxElement;
+}
+
+inline std::vector<double> vectorDiff (Vector u, Vector v)
+{
+	if (u.size() != v.size())
+		return Vector();
+
+	Vector diff (u.size(), 0);
+	for(int i=0; i<u.size(); i++)
+		diff[i] = v[i]-u[i];
+	return diff;
+}
+
+double vectorErrorL1 (Vector u, Vector v)
+{
+	auto err = vectorDiff(u,v);
+	return vectorNormL1(err);
+}
+
+double vectorErrorL2 (Vector u, Vector v)
+{
+	auto err = vectorDiff(u,v);
+	return vectorNormL2(err);
+}
+
+double vectorErrorLInf (Vector u, Vector v)
+{
+	auto err = vectorDiff(u,v);
+	return vectorNormLInf(err);
+}
+
+double dotProduct (Vector u, Vector v)
+{
+	double sum = 0;
+	for (int i=0;i<u.size();++i)
+		sum += u[i]*v[i];
+	return sum;
+}
+
+void printVector (Vector v)
+{
+	if (v.size() == 0)
+		return;
+
+	std::cout << "[" << v[0];
+	for (int i=1; i<v.size(); i++)
+		std::cout << ", " << v[i];
+	std::cout << "]" << std::endl;
+}
diff --git a/Gauss-Elim/vectorCode.h b/Gauss-Elim/vectorCode.h
new file mode 100644
index 0000000..349daf1
--- /dev/null
+++ b/Gauss-Elim/vectorCode.h
@@ -0,0 +1,17 @@
+#ifndef VECTOR_CODE_H
+#define VECTOR_CODE_H
+
+#include<vector>
+
+typedef std::vector<double> Vector;
+
+double vectorNormL1 (Vector v);
+double vectorNormL2 (Vector v);
+double vectorNormLInf (Vector v);
+double vectorErrorL1 (Vector u, Vector v);
+double vectorErrorL2 (Vector u, Vector v);
+double vectorErrorLInf (Vector u, Vector v);
+double dotProduct (Vector u, Vector v);
+void printVector (Vector v);
+
+#endif
diff --git a/linearAlgebra/gaussElim.cpp b/linearAlgebra/gaussElim.cpp
new file mode 100644
index 0000000..9102572
--- /dev/null
+++ b/linearAlgebra/gaussElim.cpp
@@ -0,0 +1,89 @@
+// GEtest.cpp
+// Author: Isaac Stone
+// Language: C++
+//
+// Synopsis: This .cpp file implements Guassin Elimination with
+// back-substitution for an augmented matrix. The function
+// gaussElim() takes a Matrix ( matrix.h ) as an argument and returns a
+// solution to Ax=b in the form of a Vector ( vectorCode.h ). If no solution or 
+// only the trivial solution exists, a message is printed 
+// to the screen and the program exits.
+
+#include <vector>
+#include <cstdio>
+#include <cstdlib>
+#include "matrix.h"
+#include "vectorCode.h"
+
+using namespace std;
+
+// Elimination step
+void Elim( double x, Vector &a, Vector &b ) {
+	
+  for ( unsigned i = 0; i < a.size(); i++ )
+	  a[i] = a[i] - ( x * b[i] );
+}
+
+Vector gaussElim( Matrix A ) {
+	
+  int i, j, p, n;
+  double sum, xn, x;
+  Vector tmp, sols;
+
+  x = 0.0;
+  p = 0;
+  n = A.size();
+
+  // Guassian Elimination
+  for ( i = 0; i < (n - 1) ; i++ )
+  {
+	p = i;
+	while( A[p][i] == 0 )
+	{
+	  p++;
+	  if ( p >  n )
+	  {
+		printf( "No solution or trivial solution\n" );
+		exit(1);
+	  }
+	}
+	
+	// Swap rows
+	if ( p != i ) {
+	  tmp = A[p];
+	  A[p] = A[i];
+	  A[i] = tmp;
+	}
+		
+	for ( j = i+1; j < n; j++ )
+	{
+	  x = A[j][i] / A[i][i];
+	  Elim( x, A[j], A[i] );
+	}
+  }
+
+  // Test for no solution
+  if ( A[ n - 1 ][ n - 1 ] == 0 )
+  {
+    printf( "No solution or trivial solution\n" );
+	exit(1);
+  }
+ 
+  // Perform back-substitution
+  sols.resize( n, 0 );
+  
+  xn = A[n-1][n] / A[n-1][n-1];
+  sols[n-1] = xn;
+  
+  for ( int i = n-1; i >= 0; i-- )
+  {
+	sum = 0;
+	for ( j = i+1; j < n+1; j++ ) {
+	  sum += ( A[i][j] * sols[j] );
+	}
+	xn = ( A[i][n] - sum ) / A[i][i];
+	sols[i] = xn;
+  }
+
+  return sols;
+}
diff --git a/linearAlgebra/gaussElim.h b/linearAlgebra/gaussElim.h
new file mode 100644
index 0000000..d7b720e
--- /dev/null
+++ b/linearAlgebra/gaussElim.h
@@ -0,0 +1,19 @@
+// GEtest.cpp
+// Author: Isaac Stone
+// Language: C++
+//
+// Synopsis: This file includes headers gaussElim.cpp to implement
+// Guassin Elimination with back-substitution for an augmented matrix. The function
+// gaussElim() takes a Matrix ( matrix.h ) as an argument and returns a
+// solution to Ax=b in the form of a Vector ( vectorCode.h ). If no solution or 
+// only the trivial solution exists, a message is printed 
+// to the screen and the program exits.
+//
+// Check "gaassElim.cpp" for more information.
+
+#include <vector>
+#include <cstdlib>
+#include "matrix.h"
+#include "vectorCode.h"
+
+Vector gaussElim( Matrix A );