diff --git a/Population Growth.ipynb b/Population Growth.ipynb
index 6ae6009..f232dbd 100644
--- a/Population Growth.ipynb	
+++ b/Population Growth.ipynb	
@@ -225,6 +225,122 @@
     "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": 7,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "from matrix_math import *"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "leslie = Matrix([[0,0.25,0.6,0.8,0.15,0],\n",
+    "[0.7,0,0,0,0,0],\n",
+    "[0,0.95,0,0,0,0],\n",
+    "[0,0,0.9,0,0,0],\n",
+    "[0,0,0,0.9,0,0],\n",
+    "[0,0,0,0,0.5,0]])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Vector([[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]])"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "leslie"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Vector([10, 0, 0, 0, 0, 0])"
+      ]
+     },
+     "execution_count": 25,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "drakes = Vector([10,0,0,0,0,0])\n",
+    "drakes"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 34,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[10, 0, 0, 0, 0, 0] Drakes in year 0 \n",
+      "[0.0, 7.0, 0.0, 0.0, 0.0, 0.0] Drakes in year 2 \n",
+      "[1.0, 0.0, 6.0, 0.0, 0.0, 0.0] Drakes in year 4 \n",
+      "[3.0, 1.0, 0.0, 5.0, 0.0, 0.0] Drakes in year 6 \n",
+      "[5.0, 2.0, 1.0, 0.0, 5.0, 0.0] Drakes in year 8 \n",
+      "[2.0, 3.0, 2.0, 1.0, 0.0, 2.0] Drakes in year 10 \n",
+      "[3.0, 1.0, 3.0, 2.0, 0.0, 0.0] Drakes in year 12 \n",
+      "[4.0, 2.0, 1.0, 3.0, 2.0, 0.0] Drakes in year 14 \n",
+      "[4.0, 3.0, 2.0, 1.0, 2.0, 1.0] Drakes in year 16 \n",
+      "[3.0, 2.0, 2.0, 1.0, 1.0, 1.0] Drakes in year 18 \n",
+      "[4.0, 2.0, 2.0, 2.0, 1.0, 0.0] Drakes in year 20 \n",
+      "[4.0, 3.0, 2.0, 2.0, 2.0, 0.0] Drakes in year 22 \n",
+      "[4.0, 3.0, 2.0, 2.0, 2.0, 1.0] Drakes in year 24 \n",
+      "[4.0, 3.0, 3.0, 2.0, 1.0, 1.0] Drakes in year 26 \n",
+      "[5.0, 3.0, 3.0, 2.0, 2.0, 0.0] Drakes in year 28 \n",
+      "[5.0, 3.0, 3.0, 2.0, 2.0, 1.0] Drakes in year 30 \n",
+      "[5.0, 3.0, 3.0, 2.0, 2.0, 1.0] Drakes in year 32 \n",
+      "[5.0, 3.0, 3.0, 3.0, 2.0, 1.0] Drakes in year 34 \n",
+      "[5.0, 3.0, 3.0, 3.0, 2.0, 1.0] Drakes in year 36 \n",
+      "[5.0, 4.0, 3.0, 3.0, 2.0, 1.0] Drakes in year 38 \n",
+      "[6.0, 4.0, 3.0, 3.0, 2.0, 1.0] Drakes in year 40 \n"
+     ]
+    }
+   ],
+   "source": [
+    "drakes = Vector([10,0,0,0,0,0])\n",
+    "multiplier = leslie\n",
+    "for i in range(0,41,2):\n",
+    "    print('{} Drakes in year {} '.format([drake//1 for drake in drakes.my_list], i))\n",
+    "    drakes = leslie * drakes\n",
+    "    \n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
diff --git a/README.md b/README.md
index 97f003e..4fbe5e3 100644
--- a/README.md
+++ b/README.md
@@ -1,3 +1,8 @@
+This module may be used with "import matrix_math" within a python script.
+
+test_matrix_math.py is nose and py.test compatible (run with 'nosetests' or 'py.test' respectively)
+
+===========================
 # Vector and Matrix Objects
 
 ## Description
diff --git a/matrix_math.py b/matrix_math.py
new file mode 100644
index 0000000..efc79d5
--- /dev/null
+++ b/matrix_math.py
@@ -0,0 +1,112 @@
+class ShapeException(Exception):
+    pass
+
+
+class Vector():
+    def __init__(self, my_list):
+        self.my_list = my_list
+        self.my_shape = self.shape()
+
+    def __add__(self, other):
+        if self.shape() != other.shape():
+            raise ShapeException
+        return Vector([self.my_list[i] + other.my_list[i]
+                       for i in range(self.shape()[0])])
+
+    def __sub__(self, other):
+        other = other * -1
+        return self + other
+
+    def __mul__(self, other):
+        if type(other) is int or type(other) is float:
+            return Vector([i * other for i in self.my_list])
+        else:
+            raise ValueError('Vector can only be multiplied by \
+                              scalar (int/float)')
+
+    def __eq__(self, other):
+        return self.my_list == other.my_list
+
+    def __repr__(self):
+        return 'Vector({})'.format(self.my_list)
+
+    def shape(self):
+        try:
+            if len(self.my_list[0]) > 1:
+                raise ValueError('Not a 1d vector')
+        except TypeError:
+            return (len(self.my_list),)
+        except:
+            raise ValueError('Not a 1d vector')
+
+    def dot(self, other):
+        if self.shape() != other.shape():
+            raise ShapeException
+        return sum([self.my_list[i] * other.my_list[i]
+                    for i in range(len(self.my_list))])
+
+    def magnitude(self):
+        return self.dot(self)**0.5
+
+
+class Matrix():
+    def __init__(self, my_list):
+        self.my_list = my_list
+        self.my_shape = self.shape()
+
+    @classmethod
+    def create(cls, size=(2, 2), fn = lambda x, y: 1):
+        return Matrix([[fn(i, j) for j in range(size[0])]
+                       for i in range(size[1])])
+
+    def __add__(self, other):
+        if self.shape() != other.shape():
+            raise ShapeException
+        return Matrix([self.my_list[i] + other.my_list[i]
+                       for i in range(self.shape()[0])])
+
+    def __sub__(self, other):
+        other = other * -1
+        return self + other
+
+    def __mul__(self, other):
+        if type(other) is int or type(other) is float:
+            return Matrix([[i * other for i in row] for row in self.my_list])
+
+        elif type(other) is Vector:
+            if self.shape()[1] != other.shape()[0]:
+                raise ShapeException
+
+            step1 = [[val * other.my_list[idx] for idx, val in enumerate(row)]
+                     for row in self.my_list]
+            return Vector([sum(x) for x in step1])
+
+        elif type(other) is Matrix:
+            if self.shape()[1] != other.shape()[0]:
+                raise ShapeException
+
+            y_transposed = [[row[i] for row in other.my_list]
+                            for i in range(len(other.my_list[0]))]
+
+            return Matrix([[Vector(row).dot(Vector(col))
+                            for col in y_transposed] for row in self.my_list])
+
+        else:
+            raise ValueError('Matrix can only be multiplied by '
+                               'scalar, Vector, or Matrix')
+
+    def __eq__(self, other):
+        return self.my_list == other.my_list
+
+    def __repr__(self):
+        return 'Matrix({})'.format(self.my_list)
+
+    def shape(self):
+        try:
+            return (len(self.my_list), len(self.my_list[0]))
+        except (TypeError, IndexError):
+            raise ValueError('Not a 2d matrix')
+
+    def rotate(self):
+        tuple_list = list(zip(*self.my_list[::-1]))
+        return Matrix([list(row) for row in tuple_list])
diff --git a/nt-coverage.sh b/nt-coverage.sh
new file mode 100755
index 0000000..c49e516
--- /dev/null
+++ b/nt-coverage.sh
@@ -0,0 +1,2 @@
+nosetests --with-coverage --cover-html --cover-package=matrix_math
+open cover/index.html 
diff --git a/test_matrix_math.py b/test_matrix_math.py
new file mode 100644
index 0000000..5b9cd39
--- /dev/null
+++ b/test_matrix_math.py
@@ -0,0 +1,278 @@
+from matrix_math import Matrix, Vector, ShapeException
+import math
+from nose.tools import raises
+import random
+
+
+def test_vector_class():
+    v = Vector([1, 5])
+    assert isinstance(v, Vector)
+
+
+def test_matrix_class():
+    m = Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
+    assert isinstance(m, Matrix)
+
+
+@raises(ValueError)
+def test_2d_vector_exception():
+    v = Vector([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
+
+
+@raises(ValueError)
+def test_empty_vector_exception():
+    v = Vector([])
+
+
+@raises(ValueError)
+def test_1d_matrix_exception():
+    m = Matrix([1, 5])
+
+
+@raises(ValueError)
+def test_empty_matrix_exception():
+    v = Matrix([])
+
+
+m = Vector([3, 4])
+n = Vector([5, 0])
+
+v = Vector([1, 3, 0])
+w = Vector([0, 2, 4])
+u = Vector([1, 1, 1])
+y = Vector([10, 20, 30])
+z = Vector([0, 0, 0])
+
+A = Matrix([[1, 0, 0],
+            [0, 1, 0],
+            [0, 0, 1]])
+B = Matrix([[1, 2, 3],
+            [4, 5, 6],
+            [7, 8, 9]])
+C = Matrix([[1, 2],
+            [2, 1],
+            [1, 2]])
+D = Matrix([[1, 2, 3],
+            [3, 2, 1]])
+
+
+def test_matrix_class_and_shape():
+    assert m.shape() == (2,)
+    assert v.shape() == (3,)
+
+
+def test_matrix_shape():
+    m = Matrix([[0, 1, 0], [1, 0, 0], [0, 0, 1]])
+    assert m.shape() == (3, 3)
+    m = Matrix([[1, 2], [3, 4], [5, 6]])
+    assert m.shape() == (3, 2)
+    assert A.shape() == (3, 3)
+    assert C.shape() == (3, 2)
+    assert D.shape() == (2, 3)
+
+
+def test_vector_eq():
+    p = Vector([3, 4])
+    assert m == p
+
+
+def test_vector_add():
+    """
+    [a b]  + [c d]  = [a+c b+d]
+
+    Matrix + Matrix = Matrix
+    """
+    assert v + w == Vector([1, 5, 4])
+    assert u + y == Vector([11, 21, 31])
+    assert u + z == u
+
+
+def test_vector_add_is_communicative():
+    assert w + y == y + w
+
+
+def test_vector_sub():
+    """
+    [a b]  - [c d]  = [a-c b-d]
+
+    Matrix + Matrix = Matrix
+    """
+    assert v - w == Vector([1, 1, -4])
+    assert w - v == Vector([-1, -1, 4])
+    assert y - z == y
+    assert w - u == (z - (u - w))
+
+
+def test_hard_mode_vector():
+    Vector([1, 2]) + Vector([0, 4])
+    Vector([1, 2]) - Vector([0, 4])
+    Vector([1, 2]) * 3
+
+    assert Vector([1, 2]) == Vector([1, 2])  # results in True
+
+
+@raises(ShapeException)
+def test_vector_sub_checks_shapes():
+    """Shape rule: the vectors must be the same size."""
+    m - v
+
+
+def test_dot():
+    """
+    dot([a b], [c d])   = a * c + b * d
+
+    dot(Vector, Vector) = Scalar
+    """
+    assert w.dot(y) == 160
+    assert m.dot(n) == 15
+    assert u.dot(z) == 0
+
+
+@raises(ShapeException)
+def test_dot_checks_shapes():
+    """Shape rule: the vectors must be the same size."""
+    v.dot(m)
+
+
+def test_vector_multiply():
+    """
+    [a b]  *  Z     = [a*Z b*Z]
+
+    Vector * Scalar = Vector
+    """
+    assert v * 0.5 == Vector([0.5, 1.5, 0])
+    assert m * 2 == Vector([6, 8])
+
+
+def test_magnitude():
+    """
+    magnitude([a b])  = sqrt(a^2 + b^2)
+
+    magnitude(Vector) = Scalar
+    """
+    assert m.magnitude() == 5
+    assert v.magnitude() == math.sqrt(10)
+    assert y.magnitude() == math.sqrt(1400)
+    assert z.magnitude() == 0
+
+
+def test_shape_matrices():
+    """shape should take a vector or matrix and return a tuple with the
+    number of rows (for a vector) or the number of rows and columns
+    (for a matrix.)"""
+    assert A.shape() == (3, 3)
+    assert C.shape() == (3, 2)
+    assert D.shape() == (2, 3)
+
+
+def test_matrix_scalar_multiply():
+    """
+    [[a b]   *  Z   =   [[a*Z b*Z]
+     [c d]]              [c*Z d*Z]]
+
+    Matrix * Scalar = Matrix
+    """
+    assert C * 3 == Matrix([[3, 6],
+                            [6, 3],
+                            [3, 6]])
+
+
+def test_matrix_vector_multiply():
+    """
+    [[a b]   *  [x   =   [a*x+b*y
+     [c d]       y]       c*x+d*y
+     [e f]                e*x+f*y]
+
+    Matrix * Vector = Vector
+    """
+    assert A * Vector([2, 5, 4]) == Vector([2, 5, 4])
+    assert B * Vector([1, 2, 3]) == Vector([14, 32, 50])
+    assert C * Vector([3, 4]) == Vector([11, 10, 11])
+    assert D * Vector([0, 1, 2]) == Vector([8, 4])
+
+
+@raises(ShapeException)
+def test_matrix_vector_multiply_checks_shapes():
+    """Shape Rule: The number of rows of the vector must equal the number of
+    columns of the matrix."""
+    C * Vector([1, 2, 3])
+
+
+def test_matrix_matrix_multiply():
+    """
+    [[a b]   *  [[w x]   =   [[a*w+b*y a*x+b*z]
+     [c d]       [y z]]       [c*w+d*y c*x+d*z]
+     [e f]                    [e*w+f*y e*x+f*z]]
+
+    Matrix * Matrix = Matrix
+    """
+    assert A * B == B
+    assert B * C == Matrix([[8, 10],
+                            [20, 25],
+                            [32, 40]])
+    assert C * D == Matrix([[7, 6, 5],
+                            [5, 6, 7],
+                            [7, 6, 5]])
+    assert D * C == Matrix([[8, 10], [8, 10]])
+
+
+@raises(ShapeException)
+def test_matrix_matrix_multiply_checks_shapes():
+    """Shape Rule: The number of columns of the first matrix must equal the
+    number of rows of the second matrix."""
+    A * D
+
+
+def test_hard_mode_matrix():
+    Matrix([[0, 1], [1, 0]]) + Matrix([[1, 1], [0, 0]])
+    Matrix([[0, 1], [1, 0]]) - Matrix([[1, 1], [0, 0]])
+    Matrix([[0, 1], [1, 0]]) * 3
+    Matrix([[0, 1], [1, 0]]) * Vector([1, 2])
+    Matrix([[1, 1, 1], [0, 0, 0]]) * Matrix([[1, 1], [2, 2], [3, 3]])
+
+    assert not Matrix([[0, 1], [1, 0]]) == Matrix([[1, 1], [0, 0]])
+    # results in False
+
+
+def test_create_class_method():
+    m1 = Matrix.create()
+    assert m1 == Matrix([[1, 1], [1, 1]])
+    m2 = Matrix.create((3, 3), lambda x, y: x + y)
+    assert m2 == Matrix([[0, 1, 2], [1, 2, 3], [2, 3, 4]])
+    random.seed(0)
+
+    def rnd(x, y):
+        return random.randint(0, 9)
+    m3 = Matrix.create((3, 2), rnd)
+    assert m3 != Matrix([[6, 6], [0, 4], [8, 7]])
+    m2 = Matrix.create((3, 3), lambda x, y: x ** 2 + y ** 2)
+    assert m2 == Matrix([[0, 1, 4], [1, 2, 5], [4, 5, 8]])
+
+
+def test_matrix_rotate():
+    assert A.rotate() == Matrix([[0, 0, 1],
+                                 [0, 1, 0],
+                                 [1, 0, 0]])
+
+
+# Add tests to get coverage > 93%
+@raises(ValueError)
+def test_vector_bad_multiplication_is_none():
+    assert m * '1' is None
+
+@raises(ValueError)
+def test_matrix_bad_multiplication_is_none():
+    assert A * '1' is None
+
+def test_vector_matrix_repr():
+    assert 'Vector(' in repr(m)
+    assert 'Matrix(' in repr(A)
+
+@raises(ShapeException)
+def test_vector_shape_exceptions_add():
+    m + v
+
+
+@raises(ShapeException)
+def test_matrix_shape_exceptions_add():
+    A + C