diff --git a/Population Growth.ipynb b/Population Growth.ipynb
index 6ae6009..209060d 100644
--- a/Population Growth.ipynb	
+++ b/Population Growth.ipynb	
@@ -225,6 +225,152 @@
     "Using the above math plus your matrix library, calculate how many drakes of each age I will have in 20 2-year cycles."
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "from matrix_math import ShapeException, Vector, Matrix"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 77,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "leslie = Matrix([[0, 0.25, 0.6, 0.8, 0.15, 0], [0.7, 0, 0, 0, 0, 0], [0, 0.95, 0, 0, 0, 0], \\\n",
+    "                [0, 0, 0.9, 0, 0, 0], [0, 0, 0, 0.9, 0, 0], [0, 0, 0, 0, 0.5, 0]])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 78,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[[0, 0.25, 0.6, 0.8, 0.15, 0], [0.7, 0, 0, 0, 0, 0], [0, 0.95, 0, 0, 0, 0], [0, 0, 0.9, 0, 0, 0], [0, 0, 0, 0.9, 0, 0], [0, 0, 0, 0, 0.5, 0]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(leslie.values)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 79,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "dragons = Vector([10, 0, 0, 0, 0, 0])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 80,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "for i in range(20):\n",
+    "    dragons = leslie * dragons"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 81,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[6.149914942919132, 4.182186803339321, 3.8477209681295763, 3.283445777563774, 2.841646230444824, 1.4082839376892242]\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(dragons.values)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 82,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "test = Matrix([[.25, .5],[0, .95]])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 83,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "leslie20 = leslie"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 84,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[[0.6407182838330322, 0.9533865528483813, 0.8694631478701617, 0.5726496552122178, 0.09224872414378701, 0.0], [0.4304940460043394, 0.6407182838330322, 0.589207447887139, 0.38925308434056616, 0.06273280205008984, 0.0], [0.3973077463172356, 0.5842419195773176, 0.5361636137495489, 0.35706936411382273, 0.05771581452194364, 0.0], [0.34629488713166195, 0.5108242452650171, 0.4623621114070214, 0.3053003556617743, 0.04925168666345661, 0.0], [0.2955101199807398, 0.445236283454994, 0.4061729376245585, 0.265355364753195, 0.042624693456672365, 0.0], [0.1420823115222412, 0.21107865712909984, 0.19694480404940706, 0.13093009099881617, 0.021124259065338366, 0.0]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "for _ in range(20):\n",
+    "    leslie20 *= leslie\n",
+    "print(leslie20.values)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 89,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[6.407182838330322, 4.304940460043394, 3.973077463172356, 3.4629488713166197, 2.9551011998073977, 1.420823115222412]\n"
+     ]
+    }
+   ],
+   "source": [
+    "print((leslie20 * Vector([10,0,0,0,0,0])).values)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
diff --git a/matrix_math.py b/matrix_math.py
new file mode 100644
index 0000000..d7d35b0
--- /dev/null
+++ b/matrix_math.py
@@ -0,0 +1,177 @@
+import math
+
+
+class ShapeException(Exception):
+    pass
+
+
+class Vector:
+
+    def __init__(self, values):
+        self.values = self.is_vector(values)
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    @property
+    def shape(self):
+        return len(self.values),
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    def __eq__(self, other):
+        return self.values == other.values
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    def __add__(self, other):
+        if self.shape != other.shape:
+            raise ShapeException()
+        else:
+            return Vector([self.values[index] + other.values[index]
+                           for index in range(len(self.values))])
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    def __sub__(self, other):
+        if self.shape != other.shape:
+            raise ShapeException()
+
+        else:
+            return Vector([self.values[index] - other.values[index]
+                           for index in range(len(self.values))])
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    def __mul__(self, other):
+        if isinstance(other, Matrix):
+            return Vector([self.dot(Vector(row)) for row in other.values])
+        else:
+            return Vector([self.values[index] * other
+                          for index in range(len(self.values))])
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    def dot(self, other):
+        """takes the dot product of the vector and another vector"""
+        if self.shape != other.shape:
+            raise ShapeException()
+        else:
+            return sum([self.values[index] * other.values[index]
+                        for index in range(len(self.values))])
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    def mag(self):
+        """Returns the magnitude of the vector"""
+        return math.sqrt(sum([coefficient ** 2
+                              for coefficient in self.values]))
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    def is_vector(self, a_list):
+        """Checks to ensure the vector is initialized with a list of numbers"""
+        if type(a_list) != list:
+            raise ValueError(
+                "Must make Vector with numerical list")
+        else:
+            for item in a_list:
+                if not isinstance(item, (int, float)):
+                    raise ValueError(
+                        "Must make Vector with numerical list")
+            return a_list
+
+
+###############################################################################
+
+
+class Matrix:
+
+    def __init__(self, values):
+        self.values = self.is_matrix(values)
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    @property
+    def shape(self):
+        return len(self.values), len(self.values[0])
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    def __eq__(self, other):
+        return self.values == other.values
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    def __add__(self, other):
+        if self.shape != other.shape:
+            raise ShapeException()
+        else:
+            return Matrix([[self.values[row][index] + other.values[row][index]
+                            for index in range(len(self.values[0]))]
+                           for row in range(len(self.values))])
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    def __sub__(self, other):
+        if self.shape != other.shape:
+            raise ShapeException()
+        else:
+            return Matrix([[self.values[row][index] - other.values[row][index]
+                            for index in range(len(self.values[0]))]
+                           for row in range(len(self.values))])
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    def __mul__(self, other):
+        """Checks if the matrix is being multiplied by a scalar, vector,
+        or other matrix and performs the corresponding operation."""
+        if isinstance(other, (int, float)):
+            return Matrix([[self.values[row][index] * other
+                            for index in range(len(self.values[0]))]
+                           for row in range(len(self.values))])
+
+        elif isinstance(other, Vector):
+            return Vector([other.dot(Vector(row)) for row in self.values])
+
+        elif isinstance(other, Matrix):
+            return Matrix([(other.transpose() * Vector(row)).values
+                           for row in self.values])
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    def transpose(self):
+        """Returns a matrix where the rows and columns of the current matrix
+        have been transposed"""
+        return Matrix([[row[index]
+                        for row in self.values]
+                       for index in range(len(self.values[0]))])
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    def is_matrix(self, a_list):
+        """Checks to ensure the matrix is initialized with a list of
+        lists of numbers"""
+        if type(a_list) != list:
+            raise ValueError(
+                "Must make Matrix w/list of numerical lists")
+        else:
+            for index in range(len(a_list)):
+                if type(a_list[index]) != list or \
+                        len(a_list[index]) != len(a_list[(index - 1)]):
+                    raise ValueError(
+                        "Must make Matrix w/list of numerical lists")
+                else:
+                    for value in a_list[index]:
+                        if not isinstance(value, (int, float)):
+                            raise ValueError(
+                                "Must make Matrix w/list of numerical lists")
+        return a_list
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
+
+    @classmethod
+    def from_function(cls, rows, columns, function):
+        """Allows creation of a matrix from the number of rows, number of
+        columns, and an identity function."""
+        return cls([[function(x, y) for y in range(columns)]
+                    for x in range(rows)])
diff --git a/test_matrix_math.py b/test_matrix_math.py
new file mode 100644
index 0000000..a8eda11
--- /dev/null
+++ b/test_matrix_math.py
@@ -0,0 +1,228 @@
+from matrix_math import ShapeException, Vector, Matrix
+from nose.tools import raises
+import math
+
+
+def test_vector_has_content():
+    v1 = Vector([1, 2, 3])
+    assert v1.values == [1, 2, 3]
+
+
+def test_vector_has_shape():
+    v1 = Vector([1, 2, 3])
+    v2 = Vector([1, 2])
+    assert v1.shape == (3,)
+    assert v2.shape == (2,)
+
+
+def test_vector_equality():
+    v1 = Vector([1, 2, 3])
+    v2 = Vector([1, 2, 3])
+    assert v1 == v2
+
+
+def test_can_add_vector_objects():
+    v1 = Vector([1, 2, 3])
+    v2 = Vector([4, 5, 6])
+    assert v1 + v2 == Vector([5, 7, 9])
+    assert v2 + v1 == Vector([5, 7, 9])
+
+
+@raises(ShapeException)
+def test_vector_add_raises_shape_exception():
+    v1 = Vector([1, 2, 3])
+    v2 = Vector([1, 2])
+    v1 + v2
+
+
+def test_can_subtract_vector_objects():
+    v1 = Vector([1, 2, 3])
+    v2 = Vector([4, 5, 6])
+    assert v1 - v2 == Vector([-3, -3, -3])
+    assert v2 - v1 == Vector([3, 3, 3])
+
+
+@raises(ShapeException)
+def test_vector_subtract_raises_shape_exception():
+    v1 = Vector([1, 2, 3])
+    v2 = Vector([1, 2])
+    v1 - v2
+
+
+def test_vector_multiplication():
+    v1 = Vector([1, 2, 3])
+    assert v1 * 2 == Vector([2, 4, 6])
+    assert v1 * 2 == Vector([2, 4, 6])
+
+
+def test_vector_dot():
+    v1 = Vector([1, 2, 3])
+    v2 = Vector([4, 5, 6])
+    assert v1.dot(v2) == 32
+
+
+@raises(ShapeException)
+def test_vector_dot_raises_shape_exception():
+    v1 = Vector([1, 2, 3])
+    v2 = Vector([1, 2])
+    v1.dot(v2)
+
+
+def test_vector_magnitude():
+    v1 = Vector([3, 4])
+    v2 = Vector([1, 3, 0])
+    v3 = Vector([10, 20, 30])
+    v4 = Vector([0, 0, 0])
+    assert v1.mag() == 5
+    assert v2.mag() == math.sqrt(10)
+    assert v3.mag() == math.sqrt(1400)
+    assert v4.mag() == 0
+
+
+def test_vector_works_with_floats():
+    v1 = Vector([1.0, 2.0, 3.0])
+    v2 = Vector([4.0, 5.0, 6.0])
+    assert v1 + v2 == Vector([5.0, 7.0, 9.0])
+    assert v2 + v1 == Vector([5.0, 7.0, 9.0])
+
+
+@raises(ValueError)
+def test_vector_creation_list_not_all_numbers():
+    v1 = Vector([1, 2, "3"])
+
+
+@raises(ValueError)
+def test_vector_creation_not_a_list():
+    v2 = Vector("hello")
+
+
+@raises(ValueError)
+def test_vector_creation2():
+    v3 = Vector(1)
+
+m1 = Matrix([[1, 0, 0],
+             [0, 1, 0],
+             [0, 0, 1]])
+
+m2 = Matrix([[1, 2, 3],
+             [4, 5, 6],
+             [7, 8, 9]])
+
+m3 = Matrix([[1, 2],
+             [2, 1],
+             [1, 2]])
+
+m4 = Matrix([[1, 2, 3],
+             [3, 2, 1]])
+
+
+def test_matrix_has_content():
+    assert m1.values == [[1, 0, 0],
+                         [0, 1, 0],
+                         [0, 0, 1]]
+
+
+def test_matrix_has_shape():
+    assert m1.shape == (3, 3)
+    assert m3.shape == (3, 2)
+    assert m4.shape == (2, 3)
+
+
+def test_matrix_equality():
+    assert m1 == Matrix([[1, 0, 0],
+                         [0, 1, 0],
+                         [0, 0, 1]])
+
+
+def test_matrix_addition():
+    assert m1 + m2 == Matrix([[2, 2, 3], [4, 6, 6], [7, 8, 10]])
+
+
+@raises(ShapeException)
+def test_matrix_add_raises_ShapeException():
+    m1 + m3
+
+
+def test_matrix_subtraction():
+    assert m2 - m1 == Matrix([[0, 2, 3], [4, 4, 6], [7, 8, 8]])
+
+
+@raises(ShapeException)
+def test_matrix_sub_raises_ShapeException():
+    m1 - m3
+
+
+def test_matrix_scalar_multiply():
+    assert m1 * 3 == Matrix([[3, 0, 0],
+                             [0, 3, 0],
+                             [0, 0, 3]])
+    assert m3 * 3 == Matrix([[3, 6],
+                             [6, 3],
+                             [3, 6]])
+v1 = Vector([2, 5, 4])
+v2 = Vector([1, 2, 3])
+v3 = Vector([3, 4])
+v4 = Vector([0, 1, 2])
+
+
+def test_matrix_vector_multiply():
+    assert m1 * v1 == Vector([2, 5, 4])
+    assert v1 * m1 == Vector([2, 5, 4])
+    assert m2 * v2 == Vector([14, 32, 50])
+    assert m3 * v3 == Vector([11, 10, 11])
+    assert m4 * v4 == Vector([8, 4])
+
+
+@raises(ShapeException)
+def test_matrix_vector_multiply_raises_ShapeException():
+    m1 * v3
+
+
+def test_matrix_transpose():
+    assert m2.transpose() == Matrix([[1, 4, 7],
+                                     [2, 5, 8],
+                                     [3, 6, 9]])
+
+
+def test_matrix_matrix_multiply():
+    assert m1 * m2 == m2
+    assert m2 * m3 == Matrix([[8, 10],
+                              [20, 25],
+                              [32, 40]])
+    assert m3 * m4 == Matrix([[7, 6, 5],
+                              [5, 6, 7],
+                              [7, 6, 5]])
+    assert m4 * m3 == Matrix([[8, 10],
+                              [8, 10]])
+    assert Matrix([[1, 2], [3, 4]]) * Matrix([[1, 2, 3], [4, 5, 6]]) == \
+        Matrix([[9, 12, 15], [19, 26, 33]])
+
+
+@raises(ShapeException)
+def test_matrix_matrix_multiply_checks_shapes():
+    m1 * m4
+
+
+@raises(ValueError)
+def test_matrix_creation_not_all_list_items_are_numbers():
+    m1 = Matrix([[3, 4, 5], [2, 5, 7], [1, 2, "3"]])
+
+
+@raises(ValueError)
+def test_matrix_creation_not_all_lists_the_same_length():
+    m2 = Matrix([[3, 4, 5], [2], [1, 2, 3]])
+
+
+@raises(ValueError)
+def test_matrix_creation_not_a_list():
+    m2 = Matrix("hello")
+
+
+def test_from_function():
+    def diagonal_ones(x, y):
+        if x == y:
+            return 1
+        else:
+            return 0
+    assert Matrix.from_function(2, 2, diagonal_ones) == Matrix([[1, 0],
+                                                                [0, 1]])