diff --git a/README.md b/README.md
index b8b91f4..2e403cd 100644
--- a/README.md
+++ b/README.md
@@ -1,3 +1,4 @@
+<<<<<<< HEAD
 # mmathlib
 
 **Compile and run**
@@ -8,6 +9,8 @@ g++ -o fractions fractions.cpp
 ./fractions
 g++ -o cnumber cnumber.cpp
 ./cnumber
+g++ -o matrix matrix.cpp           
+./matrix
 ```
 
 **Output**
@@ -40,4 +43,16 @@ a * b =  18-i
 (a * b)* = 18+i
 a / b = (6/25)+(17/25)i
 (a / b)* = (6/25)-(17/25)i
+The matrix m:
+| 7 | 5 | i |
+| 0 | 2 | 0 |
+| -i | 0 | 4 |
+The matrix m's transpose:
+| 7 | 0 | -i |
+| 5 | 2 | 0 |
+| i | 0 | 4 |
+The matrix m is not hermitian, here is the hermitian conjugate:
+| 7 | 0 | i |
+| 5 | 2 | 0 |
+| -i | 0 | 4 |
 ```
diff --git a/matrix.cpp b/matrix.cpp
new file mode 100644
index 0000000..87eb15b
--- /dev/null
+++ b/matrix.cpp
@@ -0,0 +1,50 @@
+#include "matrix.hpp"
+#include "../complex-numbers/cnumber.hpp"
+#include "../vectors/vector.hpp"
+
+int main() {
+  vector a = vector(2);
+  a[0] = cnumber(7, 0);
+  a[1] = cnumber(6, 0);
+
+  vector b = vector(2);
+  b[0] = cnumber(5, 0);
+  b[1] = cnumber(3, 0);
+
+  vector c = vector(2);
+  c[0] = cnumber(2, 0);
+  c[1] = cnumber(5, 0);
+
+
+  vector d = vector(2);
+  d[0] = cnumber(1, 0);
+  d[1] = cnumber(0, 1);
+
+  matrix m = matrix(2, 2);
+  m[0] = a;
+  m[1] = b;
+
+  matrix n = matrix(2, 2);
+  n[0] = c;
+  n[1] = d;
+
+  cout << "The matrix m:" << endl;
+  cout << m << endl;
+  cout << "The matrix m's transpose:" << endl;
+  cout << m.transpose() << endl;
+  cout << "The matrix n" << endl;
+  cout << n << endl;
+  cout << "The matrix m * n" << endl;
+  cout << m * n << endl;
+  if  (!m.is_hermitian()) {
+    cout << "The matrix m is not hermitian, here is the hermitian conjugate:"
+         << endl;
+    cout << m.hermitian_conjugate() << endl;
+  }
+  if  (!n.is_hermitian()) {
+    cout << "The matrix n is not hermitian, here is the hermitian conjugate:"
+         << endl;
+    cout << n.hermitian_conjugate() << endl;
+  }
+  return 0;
+}
diff --git a/matrix.hpp b/matrix.hpp
new file mode 100644
index 0000000..659ee1c
--- /dev/null
+++ b/matrix.hpp
@@ -0,0 +1,171 @@
+#include "../vectors/vector.hpp"
+#include <iostream>
+#include <sstream>
+#include <string>
+using namespace std;
+class matrix {
+private:
+  long long num_cols;
+  long long num_rows;
+  vector *cols;
+
+public:
+  matrix(const long long num_cols, const long long num_rows) {
+    this->num_cols = num_rows;
+    this->num_rows = num_cols;
+    this->cols = (vector *)calloc(sizeof(vector), num_rows);
+    for (long long i = 0; i < num_rows; i++) {
+      this->cols[i] = vector(num_cols);
+    }
+  }
+  ~matrix() {
+    this->num_cols = 0;
+    this->num_rows = 0;
+    //    free(this->cols);
+    // this->cols = NULL;
+  }
+  const vector get_diagonal() const {
+    long long diag_len = this->num_rows;
+    if (this->num_cols < diag_len) {
+      diag_len = this->num_cols;
+    }
+    vector v(diag_len);
+    for (long long i = 0; i < this->num_cols; i++) {
+      for (long long j = 0; j < this->num_rows; j++) {
+        if (i == j) {
+          v[i] = this->cols[i][j];
+        }
+      }
+    }
+    return v;
+  }
+  const vector get_row(long long index) const { return this->cols[index]; }
+
+  const vector get_column(long long index) const {
+    vector v(this->num_rows);
+    for (long long j = 0; j < this->num_rows; j++) {
+      v[j] = this->cols[j][index];
+    }
+    return v;
+  }
+
+  const matrix rotate_ninety() const {
+    matrix m(this->num_cols, this->num_rows);
+    return m;
+  }
+  const matrix transpose() const {
+    matrix n(this->num_rows, this->num_cols);
+    for (long long i = 0; i < this->num_cols; i++) {
+      for (long long j = 0; j < this->num_rows; j++) {
+        n[j][i] = this->cols[i][j];
+      }
+    }
+
+    return n;
+  }
+  const matrix conjugate() const {
+    matrix n(this->num_cols, this->num_rows);
+    for (long long i = 0; i < this->num_cols; i++) {
+      for (long long j = 0; j < this->num_rows; j++) {
+        n[i][j] = this->cols[i][j].conjugate();
+      }
+    }
+
+    return n;
+  }
+  const matrix hermitian_conjugate() const {
+    return this->transpose().conjugate();
+  }
+  const bool is_hermitian() const {
+    if (this->num_rows != this->num_cols)
+      return false;
+    matrix m = this->hermitian_conjugate();
+    bool equal = true;
+    for (long long i = 0; i < m.num_cols; i++) {
+      for (long long j = 0; j < m.num_rows; j++) {
+        if (m[i][j] != this->cols[i][j]) {
+          equal = false;
+        }
+      }
+    }
+    return equal;
+  }
+  friend ostream &operator<<(ostream &os, const matrix &m) {
+    char last = '\0';
+    for (long long i = 0; i < m.num_cols; i++) {
+      for (long long j = 0; j < m.num_rows; j++) {
+        string symbols[3];
+        symbols[0] = "|";
+        ostringstream oss;
+        oss << " " << m.cols[i][j] << " ";
+        symbols[1] = oss.str();
+        symbols[2] = "|";
+        for (int i = 0; i < 3; i++) {
+          int len = symbols[i].length() - 1;
+          char cur = symbols[i][0];
+          if (cur != last) {
+            os << symbols[i];
+          }
+          last = symbols[i][len];
+        }
+      }
+      if (i != m.num_cols - 1)
+        os << endl << "|";
+    }
+    return os;
+  }
+  const bool operator==(const matrix &m) const {
+    bool equal = true;
+    for (long long i = 0; i < this->num_cols; i++) {
+      for (long long j = 0; j < this->num_rows; j++) {
+        if (this->cols[i][j] != m[i][j]) {
+          equal = false;
+        }
+      }
+    }
+    return equal;
+  }
+  const vector operator[](const long long index) const {
+    return this->cols[index];
+  }
+  const matrix operator*(const cnumber z) const {
+    matrix n(this->num_cols, this->num_rows);
+    for (long long i = 0; i < this->num_cols; i++) {
+      n[i] = this->cols[i] * z;
+    }
+
+    return n;
+  }
+  const matrix operator*(const matrix m) const {
+    matrix n(this->num_cols, m.num_rows);
+    // if (this->num_cols != m.num_rows && m.num_cols != this->num_rows)
+    //   return n;
+    for (long long i = 0; i < this->num_rows; i++) {
+      for (long long j = 0; j < this->num_rows; j++) {
+        n[i][j] = this->get_row(i) * m.get_column(j);
+      }
+    }
+
+    return n;
+  }
+
+  const matrix operator+(const matrix &m) const {
+    matrix n(this->num_cols, this->num_rows);
+    for (long long i = 0; i < this->num_cols; i++) {
+      n[i] = this->cols[i] + m[i];
+    }
+
+    return n;
+  }
+  const matrix operator-(const matrix &m) const {
+    matrix n(this->num_cols, this->num_rows);
+    for (long long i = 0; i < this->num_cols; i++) {
+      n[i] = this->cols[i] - m[i];
+    }
+
+    return n;
+  }
+  // FIXME: Figure out how to make sure you do not try to access something
+  // outside of the index
+  vector &operator[](const long long index) { return this->cols[index]; }
+};