diff --git a/doc/software_manual.org b/doc/software_manual.org
index 981d863..eebcf17 100644
--- a/doc/software_manual.org
+++ b/doc/software_manual.org
@@ -5,7 +5,7 @@
 #+STARTUP: entitiespretty fold inlineimages
 
 * Design
-The LIZFCM static library (at [[https://github.com/Simponic/math-4610][[https://github.com/Simponic/math-4610]] is a successor to my
+The LIZFCM static library (at [[https://github.com/Simponic/math-4610]]) is a successor to my
 attempt at writing codes for the Fundamentals of Computational Mathematics course in Common
 Lisp, but the effort required to meet the requirement of creating a static library became
 too difficult to integrate outside of the ~ASDF~ solution that Common Lisp already brings
@@ -14,13 +14,14 @@ to the table.
 All of the work established in ~deprecated-cl~ has been painstakingly translated into
 the C programming language. I have a couple tenets for its design:
 
-+ Implemntations of routines should all be done immutably in respect to arguments.
++ Implementations of routines should all be done immutably in respect to arguments.
 + Functional programming is good (it's... rough in C though).
-+ Routines are separated into "module" c files, and not individual files per function.
++ Routines are separated into "modules" that follow a form of separation of concerns
+  in files, and not individual files per function.
 
 * Compilation
-A provided ~Makefile~ is added for convencience. It has been tested on an M1 machine running MacOS as
-well as Arch Linux.
+A provided ~Makefile~ is added for convencience. It has been tested on an ~arm~-based M1 machine running
+MacOS as well as ~x86~ Arch Linux.
 
 1. ~cd~ into the root of the repo
 2. ~make~
@@ -317,6 +318,39 @@ void free_vector(Array_double *v) {
 }
 #+END_SRC
 
+*** ~add_element~
++ Author: Elizabeth Hunt
++ Name: ~add_element~
++ Location: ~src/vector.c~
++ Input: a pointer to an ~Array_double~
++ Output: a new ~Array_double~ with element ~x~ appended.
+
+#+BEGIN_SRC c
+Array_double *add_element(Array_double *v, double x) {
+  Array_double *pushed = InitArrayWithSize(double, v->size + 1, 0.0);
+  for (size_t i = 0; i < v->size; ++i)
+    pushed->data[i] = v->data[i];
+  pushed->data[v->size] = x;
+  return pushed;
+}
+#+END_SRC
+
+*** ~slice_element~
++ Author: Elizabeth Hunt
++ Name: ~slice_element~
++ Location: ~src/vector.c~
++ Input: a pointer to an ~Array_double~
++ Output: a new ~Array_double~ with element ~x~ sliced.
+
+#+BEGIN_SRC c
+Array_double *slice_element(Array_double *v, size_t x) {
+  Array_double *sliced = InitArrayWithSize(double, v->size - 1, 0.0);
+  for (size_t i = 0; i < v->size - 1; ++i)
+    sliced->data[i] = i >= x ? v->data[i + 1] : v->data[i];
+  return sliced;
+}
+#+END_SRC
+
 *** ~copy_vector~
 + Author: Elizabeth Hunt
 + Name: ~copy_vector~
@@ -370,55 +404,54 @@ void format_vector_into(Array_double *v, char *s) {
   matrix $L$, $U$, respectively such that $LU = m$.
 + Output: a pointer to the location in memory in which two ~Matrix_double~'s reside: the first
   representing $L$, the second, $U$.
-+ Errors: Exits and throws a status code of ~-1~ when encountering a matrix that cannot be
++ Errors: Fails assertions when encountering a matrix that cannot be
   decomposed
 
 #+BEGIN_SRC c
-  Matrix_double **lu_decomp(Matrix_double *m) {
-    assert(m->cols == m->rows);
+Matrix_double **lu_decomp(Matrix_double *m) {
+  assert(m->cols == m->rows);
 
-    Matrix_double *u = copy_matrix(m);
-    Matrix_double *l_empt = InitMatrixWithSize(double, m->rows, m->cols, 0.0);
-    Matrix_double *l = put_identity_diagonal(l_empt);
-    free(l_empt);
+  Matrix_double *u = copy_matrix(m);
+  Matrix_double *l_empt = InitMatrixWithSize(double, m->rows, m->cols, 0.0);
+  Matrix_double *l = put_identity_diagonal(l_empt);
+  free_matrix(l_empt);
 
+  Matrix_double **u_l = malloc(sizeof(Matrix_double *) * 2);
 
-    Matrix_double **u_l = malloc(sizeof(Matrix_double *) * 2);
-
-    for (size_t y = 0; y < m->rows; y++) {
-      if (u->data[y]->data[y] == 0) {
-        printf("ERROR: a pivot is zero in given matrix\n");
-        exit(-1);
-      }
+  for (size_t y = 0; y < m->rows; y++) {
+    if (u->data[y]->data[y] == 0) {
+      printf("ERROR: a pivot is zero in given matrix\n");
+      assert(false);
     }
-
-    if (u && l) {
-      for (size_t x = 0; x < m->cols; x++) {
-        for (size_t y = x + 1; y < m->rows; y++) {
-          double denom = u->data[x]->data[x];
-
-          if (denom == 0) {
-            printf("ERROR: non-factorable matrix\n");
-            exit(-1);
-          }
-
-          double factor = -(u->data[y]->data[x] / denom);
-
-          Array_double *scaled = scale_v(u->data[x], factor);
-          Array_double *added = add_v(scaled, u->data[y]);
-          free_vector(scaled);
-          free_vector(u->data[y]);
-
-          u->data[y] = added;
-          l->data[y]->data[x] = -factor;
-        }
-      }
-    }
-
-    u_l[0] = u;
-    u_l[1] = l;
-    return u_l;
   }
+
+  if (u && l) {
+    for (size_t x = 0; x < m->cols; x++) {
+      for (size_t y = x + 1; y < m->rows; y++) {
+        double denom = u->data[x]->data[x];
+
+        if (denom == 0) {
+          printf("ERROR: non-factorable matrix\n");
+          assert(false);
+        }
+
+        double factor = -(u->data[y]->data[x] / denom);
+
+        Array_double *scaled = scale_v(u->data[x], factor);
+        Array_double *added = add_v(scaled, u->data[y]);
+        free_vector(scaled);
+        free_vector(u->data[y]);
+
+        u->data[y] = added;
+        l->data[y]->data[x] = -factor;
+      }
+    }
+  }
+
+  u_l[0] = u;
+  u_l[1] = l;
+  return u_l;
+}
 #+END_SRC
 *** ~bsubst~
 + Author: Elizabeth Hunt
@@ -467,7 +500,7 @@ Array_double *fsubst(Matrix_double *l, Array_double *b) {
 }
 #+END_SRC
 
-*** ~solve_matrix~
+*** ~solve_matrix_lu_bsubst~
 + Author: Elizabeth Hunt
 + Location: ~src/matrix.c~
 + Input: a pointer to a ~Matrix_double~ $m$ and a pointer to an ~Array_double~ $b$
@@ -480,7 +513,7 @@ Here we make use of forward substitution to first solve $Ly = b$ given $L$ as th
 Then, $LUx = b$, thus $x$ is a solution.
 
 #+BEGIN_SRC c
-Array_double *solve_matrix(Matrix_double *m, Array_double *b) {
+Array_double *solve_matrix_lu_bsubst(Matrix_double *m, Array_double *b) {
   assert(b->size == m->rows);
   assert(m->rows == m->cols);
 
@@ -495,11 +528,97 @@ Array_double *solve_matrix(Matrix_double *m, Array_double *b) {
 
   free_matrix(u);
   free_matrix(l);
+  free(u_l);
 
   return x;
 }
 #+END_SRC
 
+*** ~gaussian_elimination~
++ Author: Elizabeth Hunt
++ Location: ~src/matrix.c~
++ Input: a pointer to a ~Matrix_double~ $m$ 
++ Output: a pointer to a copy of $m$ in reduced echelon form
+
+This works by finding the row with a maximum value in the column $k$. Then, it uses that as a pivot, and
+applying reduction to all other rows. The general idea is available at [[https://en.wikipedia.org/wiki/Gaussian_elimination]].
+
+#+BEGIN_SRC c
+Matrix_double *gaussian_elimination(Matrix_double *m) {
+  uint64_t h = 0;
+  uint64_t k = 0;
+
+  Matrix_double *m_cp = copy_matrix(m);
+
+  while (h < m_cp->rows && k < m_cp->cols) {
+    uint64_t max_row = 0;
+    double total_max = 0.0;
+
+    for (uint64_t row = h; row < m_cp->rows; row++) {
+      double this_max = c_max(fabs(m_cp->data[row]->data[k]), total_max);
+      if (c_max(this_max, total_max) == this_max) {
+        max_row = row;
+      }
+    }
+
+    if (max_row == 0) {
+      k++;
+      continue;
+    }
+
+    Array_double *swp = m_cp->data[max_row];
+    m_cp->data[max_row] = m_cp->data[h];
+    m_cp->data[h] = swp;
+
+    for (uint64_t row = h + 1; row < m_cp->rows; row++) {
+      double factor = m_cp->data[row]->data[k] / m_cp->data[h]->data[k];
+      m_cp->data[row]->data[k] = 0.0;
+
+      for (uint64_t col = k + 1; col < m_cp->cols; col++) {
+        m_cp->data[row]->data[col] -= m_cp->data[h]->data[col] * factor;
+      }
+    }
+
+    h++;
+    k++;
+  }
+
+  return m_cp;
+}
+#+END_SRC
+
+*** ~solve_matrix_gaussian~
++ Author: Elizabeth Hunt
++ Location: ~src/matrix.c~
++ Input: a pointer to a ~Matrix_double~ $m$ and a target ~Array_double~ $b$
++ Output: a pointer to a vector $x$ being the solution to the equation $mx = b$
+
+We first perform ~gaussian_elimination~ after augmenting $m$ and $b$. Then, as $m$ is in reduced echelon form, it's an upper
+triangular matrix, so we can perform back substitution to compute $x$.
+
+#+BEGIN_SRC c
+Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b) {
+  assert(b->size == m->rows);
+  assert(m->rows == m->cols);
+
+  Matrix_double *m_augment_b = add_column(m, b);
+  Matrix_double *eliminated = gaussian_elimination(m_augment_b);
+
+  Array_double *b_gauss = col_v(eliminated, m->cols);
+  Matrix_double *u = slice_column(eliminated, m->rows);
+
+  Array_double *solution = bsubst(u, b_gauss);
+
+  free_matrix(m_augment_b);
+  free_matrix(eliminated);
+  free_matrix(u);
+  free_vector(b_gauss);
+
+  return solution;
+}
+#+END_SRC
+
+
 *** ~m_dot_v~
 + Author: Elizabeth Hunt
 + Location: ~src/matrix.c~
@@ -535,6 +654,48 @@ Matrix_double *put_identity_diagonal(Matrix_double *m) {
 }
 #+END_SRC
 
+*** ~slice_column~
++ Author: Elizabeth Hunt
++ Location: ~src/matrix.c~
++ Input: a pointer to a ~Matrix_double~
++ Output: a pointer to a copy of the given ~Matrix_double~ with column at ~x~ sliced
+
+#+BEGIN_SRC c
+Matrix_double *slice_column(Matrix_double *m, size_t x) {
+  Matrix_double *sliced = copy_matrix(m);
+
+  for (size_t row = 0; row < m->rows; row++) {
+    Array_double *old_row = sliced->data[row];
+    sliced->data[row] = slice_element(old_row, x);
+    free_vector(old_row);
+  }
+  sliced->cols--;
+
+  return sliced;
+}
+#+END_SRC
+
+*** ~add_column~
++ Author: Elizabet Hunt
++ Location: ~src/matrix.c~
++ Input: a pointer to a ~Matrix_double~ and a new vector representing the appended column ~x~
++ Output: a pointer to a copy of the given ~Matrix_double~ with a new column ~x~
+
+#+BEGIN_SRC c
+Matrix_double *add_column(Matrix_double *m, Array_double *v) {
+  Matrix_double *pushed = copy_matrix(m);
+
+  for (size_t row = 0; row < m->rows; row++) {
+    Array_double *old_row = pushed->data[row];
+    pushed->data[row] = add_element(old_row, v->data[row]);
+    free_vector(old_row);
+  }
+
+  pushed->cols++;
+  return pushed;
+}
+#+END_SRC
+
 *** ~copy_matrix~
 + Author: Elizabeth Hunt
 + Location: ~src/matrix.c~
diff --git a/doc/software_manual.pdf b/doc/software_manual.pdf
index ba0b14b..8cd70f4 100644
Binary files a/doc/software_manual.pdf and b/doc/software_manual.pdf differ
diff --git a/doc/software_manual.tex b/doc/software_manual.tex
index 292f868..c6eb910 100644
--- a/doc/software_manual.tex
+++ b/doc/software_manual.tex
@@ -1,4 +1,4 @@
-% Created 2023-10-18 Wed 13:06
+% Created 2023-10-30 Mon 09:44
 % Intended LaTeX compiler: pdflatex
 \documentclass[11pt]{article}
 \usepackage[utf8]{inputenc}
@@ -31,8 +31,8 @@
 \setlength\parindent{0pt}
 
 \section{Design}
-\label{sec:org7204e7c}
-The LIZFCM static library (at \href{https://github.com/Simponic/math-4610}{[https://github.com/Simponic/math-4610} is a successor to my
+\label{sec:org7f9c526}
+The LIZFCM static library (at \url{https://github.com/Simponic/math-4610}) is a successor to my
 attempt at writing codes for the Fundamentals of Computational Mathematics course in Common
 Lisp, but the effort required to meet the requirement of creating a static library became
 too difficult to integrate outside of the \texttt{ASDF} solution that Common Lisp already brings
@@ -42,13 +42,13 @@ All of the work established in \texttt{deprecated-cl} has been painstakingly tra
 the C programming language. I have a couple tenets for its design:
 
 \begin{itemize}
-\item Implemntations of routines should all be done immutably in respect to arguments.
+\item Implementations of routines should all be done immutably in respect to arguments.
 \item Functional programming is good (it's\ldots{} rough in C though).
 \item Routines are separated into "module" c files, and not individual files per function.
 \end{itemize}
 
 \section{Compilation}
-\label{sec:org16cc307}
+\label{sec:org911a41e}
 A provided \texttt{Makefile} is added for convencience. It has been tested on an M1 machine running MacOS as
 well as Arch Linux.
 
@@ -74,11 +74,11 @@ Which is then bundled into a static library in \texttt{lib/lizfcm.a} which can b
 in the standard method.
 
 \section{The LIZFCM API}
-\label{sec:org832532a}
+\label{sec:orgd74cd2d}
 \subsection{Simple Routines}
-\label{sec:org540b602}
+\label{sec:org66bba13}
 \subsubsection{\texttt{smaceps}}
-\label{sec:org4d03b6e}
+\label{sec:orgeae9531}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{smaceps}
@@ -104,7 +104,7 @@ float smaceps() {
 \end{verbatim}
 
 \subsubsection{\texttt{dmaceps}}
-\label{sec:org2603bfc}
+\label{sec:org237c904}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{dmaceps}
@@ -130,9 +130,9 @@ double dmaceps() {
 \end{verbatim}
 
 \subsection{Derivative Routines}
-\label{sec:org95c28e9}
+\label{sec:org9cf9027}
 \subsubsection{\texttt{central\_derivative\_at}}
-\label{sec:org950de62}
+\label{sec:org1fcd333}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{central\_derivative\_at}
@@ -163,7 +163,7 @@ double central_derivative_at(double (*f)(double), double a, double h) {
 \end{verbatim}
 
 \subsubsection{\texttt{forward\_derivative\_at}}
-\label{sec:org832eda6}
+\label{sec:org6a768fc}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{forward\_derivative\_at}
@@ -194,7 +194,7 @@ double forward_derivative_at(double (*f)(double), double a, double h) {
 \end{verbatim}
 
 \subsubsection{\texttt{backward\_derivative\_at}}
-\label{sec:org591836d}
+\label{sec:org610ce76}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{backward\_derivative\_at}
@@ -225,9 +225,9 @@ double backward_derivative_at(double (*f)(double), double a, double h) {
 \end{verbatim}
 
 \subsection{Vector Routines}
-\label{sec:org5254fe4}
+\label{sec:orgfd176e6}
 \subsubsection{Vector Arithmetic: \texttt{add\_v, minus\_v}}
-\label{sec:orgf802d61}
+\label{sec:org2dbc55f}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name(s): \texttt{add\_v}, \texttt{minus\_v}
@@ -258,7 +258,7 @@ Array_double *minus_v(Array_double *v1, Array_double *v2) {
 \end{verbatim}
 
 \subsubsection{Norms: \texttt{l1\_norm}, \texttt{l2\_norm}, \texttt{linf\_norm}}
-\label{sec:orgc56e22d}
+\label{sec:org53a2ffc}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name(s): \texttt{l1\_norm}, \texttt{l2\_norm}, \texttt{linf\_norm}
@@ -292,7 +292,7 @@ double linf_norm(Array_double *v) {
 \end{verbatim}
 
 \subsubsection{\texttt{vector\_distance}}
-\label{sec:orgb54922f}
+\label{sec:org1fb3f8c}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{vector\_distance}
@@ -313,7 +313,7 @@ double vector_distance(Array_double *v1, Array_double *v2,
 \end{verbatim}
 
 \subsubsection{Distances: \texttt{l1\_distance}, \texttt{l2\_distance}, \texttt{linf\_distance}}
-\label{sec:orgf22f8e0}
+\label{sec:org4a25a94}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name(s): \texttt{l1\_distance}, \texttt{l2\_distance}, \texttt{linf\_distance}
@@ -339,7 +339,7 @@ double linf_distance(Array_double *v1, Array_double *v2) {
 \end{verbatim}
 
 \subsubsection{\texttt{sum\_v}}
-\label{sec:org4593341}
+\label{sec:org035a547}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{sum\_v}
@@ -359,7 +359,7 @@ double sum_v(Array_double *v) {
 
 
 \subsubsection{\texttt{scale\_v}}
-\label{sec:org3123f61}
+\label{sec:org12b0853}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{scale\_v}
@@ -378,7 +378,7 @@ Array_double *scale_v(Array_double *v, double m) {
 \end{verbatim}
 
 \subsubsection{\texttt{free\_vector}}
-\label{sec:org983efcf}
+\label{sec:org70ba90c}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{free\_vector}
@@ -396,7 +396,7 @@ void free_vector(Array_double *v) {
 \end{verbatim}
 
 \subsubsection{\texttt{copy\_vector}}
-\label{sec:orgde05d32}
+\label{sec:org57afc74}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{copy\_vector}
@@ -416,7 +416,7 @@ Array_double *copy_vector(Array_double *v) {
 \end{verbatim}
 
 \subsubsection{\texttt{format\_vector\_into}}
-\label{sec:org2e779f3}
+\label{sec:orgc346c3c}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{format\_vector\_into}
@@ -446,9 +446,9 @@ void format_vector_into(Array_double *v, char *s) {
 \end{verbatim}
 
 \subsection{Matrix Routines}
-\label{sec:org2354147}
+\label{sec:org3b053ab}
 \subsubsection{\texttt{lu\_decomp}}
-\label{sec:org3690faa}
+\label{sec:org5553968}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{lu\_decomp}
@@ -509,7 +509,7 @@ Matrix_double **lu_decomp(Matrix_double *m) {
 }
 \end{verbatim}
 \subsubsection{\texttt{bsubst}}
-\label{sec:orgdeba296}
+\label{sec:org253efdc}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{bsubst}
@@ -534,7 +534,7 @@ Array_double *bsubst(Matrix_double *u, Array_double *b) {
 }
 \end{verbatim}
 \subsubsection{\texttt{fsubst}}
-\label{sec:org60d3435}
+\label{sec:orge0c7bc6}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{fsubst}
@@ -562,7 +562,7 @@ Array_double *fsubst(Matrix_double *l, Array_double *b) {
 \end{verbatim}
 
 \subsubsection{\texttt{solve\_matrix}}
-\label{sec:org914121f}
+\label{sec:orgbcd445a}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Location: \texttt{src/matrix.c}
@@ -598,7 +598,7 @@ Array_double *solve_matrix(Matrix_double *m, Array_double *b) {
 \end{verbatim}
 
 \subsubsection{\texttt{m\_dot\_v}}
-\label{sec:orgae0f4c9}
+\label{sec:orga9b1f68}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Location: \texttt{src/matrix.c}
@@ -620,7 +620,7 @@ Array_double *m_dot_v(Matrix_double *m, Array_double *v) {
 \end{verbatim}
 
 \subsubsection{\texttt{put\_identity\_diagonal}}
-\label{sec:org6d84f6a}
+\label{sec:org33ead5e}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Location: \texttt{src/matrix.c}
@@ -639,7 +639,7 @@ Matrix_double *put_identity_diagonal(Matrix_double *m) {
 \end{verbatim}
 
 \subsubsection{\texttt{copy\_matrix}}
-\label{sec:orge750c56}
+\label{sec:org34b3f5b}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Location: \texttt{src/matrix.c}
@@ -659,7 +659,7 @@ Matrix_double *copy_matrix(Matrix_double *m) {
 \end{verbatim}
 
 \subsubsection{\texttt{free\_matrix}}
-\label{sec:org4ebcf85}
+\label{sec:org9c91101}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Location: \texttt{src/matrix.c}
@@ -678,7 +678,7 @@ void free_matrix(Matrix_double *m) {
 \end{verbatim}
 
 \subsubsection{\texttt{format\_matrix\_into}}
-\label{sec:org308ee0d}
+\label{sec:org51f3e27}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{format\_matrix\_into}
@@ -705,9 +705,9 @@ void format_matrix_into(Matrix_double *m, char *s) {
 }
 \end{verbatim}
 \subsection{Root Finding Methods}
-\label{sec:org8981156}
+\label{sec:org0e83d47}
 \subsubsection{\texttt{find\_ivt\_range}}
-\label{sec:orga5835b0}
+\label{sec:org3e4e34e}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{find\_ivt\_range}
@@ -737,7 +737,7 @@ double *find_ivt_range(double (*f)(double), double start_x, double delta,
 }
 \end{verbatim}
 \subsubsection{\texttt{bisect\_find\_root}}
-\label{sec:orgb118fc7}
+\label{sec:org48f0967}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name(s): \texttt{bisect\_find\_root}
@@ -765,7 +765,7 @@ double bisect_find_root(double (*f)(double), double a, double b,
 }
 \end{verbatim}
 \subsubsection{\texttt{bisect\_find\_root\_with\_error\_assumption}}
-\label{sec:orgf5124e7}
+\label{sec:org15e3c2d}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{bisect\_find\_root\_with\_error\_assumption}
@@ -789,9 +789,9 @@ double bisect_find_root_with_error_assumption(double (*f)(double), double a,
 \end{verbatim}
 
 \subsection{Linear Routines}
-\label{sec:org6f4fce5}
+\label{sec:org98cb54b}
 \subsubsection{\texttt{least\_squares\_lin\_reg}}
-\label{sec:orge810f5f}
+\label{sec:org0c0c5d7}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Name: \texttt{least\_squares\_lin\_reg}
@@ -821,12 +821,12 @@ Line *least_squares_lin_reg(Array_double *x, Array_double *y) {
 }
 \end{verbatim}
 \subsection{Appendix / Miscellaneous}
-\label{sec:org85d2eae}
+\label{sec:orge34af18}
 \subsubsection{Data Types}
-\label{sec:org198ca2d}
+\label{sec:org0f2f877}
 \begin{enumerate}
 \item \texttt{Line}
-\label{sec:org1866885}
+\label{sec:org3f27166}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Location: \texttt{inc/types.h}
@@ -839,7 +839,7 @@ typedef struct Line {
 } Line;
 \end{verbatim}
 \item The \texttt{Array\_<type>} and \texttt{Matrix\_<type>}
-\label{sec:org4a1c956}
+\label{sec:org83fc1f3}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Location: \texttt{inc/types.h}
@@ -871,10 +871,10 @@ typedef struct {
 \end{enumerate}
 
 \subsubsection{Macros}
-\label{sec:org1976330}
+\label{sec:org2bf9bf0}
 \begin{enumerate}
 \item \texttt{c\_max} and \texttt{c\_min}
-\label{sec:org208b148}
+\label{sec:orgcaa569e}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Location: \texttt{inc/macros.h}
@@ -888,7 +888,7 @@ typedef struct {
 \end{verbatim}
 
 \item \texttt{InitArray}
-\label{sec:orgccc4528}
+\label{sec:org5805999}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Location: \texttt{inc/macros.h}
@@ -909,7 +909,7 @@ typedef struct {
 \end{verbatim}
 
 \item \texttt{InitArrayWithSize}
-\label{sec:org7e87550}
+\label{sec:org264d6b7}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Location: \texttt{inc/macros.h}
@@ -930,7 +930,7 @@ typedef struct {
 \end{verbatim}
 
 \item \texttt{InitMatrixWithSize}
-\label{sec:orge6ec2b1}
+\label{sec:org310d41a}
 \begin{itemize}
 \item Author: Elizabeth Hunt
 \item Location: \texttt{inc/macros.h}
diff --git a/homeworks/hw-5.org b/homeworks/hw-5.org
new file mode 100644
index 0000000..d650375
--- /dev/null
+++ b/homeworks/hw-5.org
@@ -0,0 +1,59 @@
+#+TITLE: Homework 5
+#+AUTHOR: Elizabeth Hunt
+#+LATEX_HEADER: \notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+#+LATEX: \setlength\parindent{0pt}
+#+OPTIONS: toc:nil
+
+* Question One
+See LIZFCM \rightarrow Matrix Routines \rightarrow ~lu decomp~ & ~bsubst~.
+
+The test ~UTEST(matrix, lu_decomp)~ is a unit test for the ~lu_decomp~ routine,
+and ~UTEST(matrix, bsubst)~ verifies back substitution on an upper triangular
+3 \times 3 matrix with a known solution that can be verified manually.
+
+Both can be found in ~tests/matrix.t.c~.
+
+* Question Two
+Unless the following are met, the resulting solution will be garbage.
+
+1. The matrix $U$ must be not be singular.
+2. $U$ must be square (or it will fail the ~assert~).
+3. The system created by $Ux = b$ must be consistent.
+4. $U$ is in (obviously) upper-triangular form.
+
+Thus, the actual calculation performing the $LU$ decomposition
+(in ~lu_decomp~) does a sanity
+check for 1-3 will fail an assert, should a point along the diagonal (pivot) be
+zero, or the matrix be non-factorable.
+
+* Question Three
+See LIZFCM \rightarrow Matrix Routines \rightarrow ~fsubst~.
+
+~UTEST(matrix, fsubst)~ verifies forward substitution on a lower triangular 3 \times 3
+matrix with a known solution that can be verified manually.
+
+* Question Four
+
+See LIZFCM \rightarrow Matrix Routines \rightarrow ~gaussian_elimination~ and ~solve_gaussian_elimination~.
+
+* Question Five
+See LIZFCM \rightarrow Matrix Routines \rightarrow ~m_dot_v~, and the ~UTEST(matrix, m_dot_v)~ in
+~tests/matrix.t.c~.
+
+* Question Six
+See ~UTEST(matrix, solve_gaussian_elimination)~ in ~tests/matrix.t.c~, which generates a diagonally dominant 10 \times 10 matrix
+and shows that the solution is consistent with the initial matrix, according to the steps given. Then,
+we do a dot product between each row of the diagonally dominant matrix and the solution vector to ensure
+it is near equivalent to the input vector.
+
+* Question Seven
+See ~UTEST(matrix, solve_matrix_lu_bsubst)~ which does the same test in Question Six with the solution according to
+~solve_matrix_lu_bsubst~ as shown in the Software Manual.
+
+* Question Eight
+No, since the time complexity for Gaussian Elimination is always less than that of the LU factorization solution by $O(n^2)$ operations
+(in LU factorization we perform both backwards and forwards substitutions proceeding the LU decomp, in Gaussian Elimination we only need
+back substitution).
+
+* Question Nine, Ten
+See LIZFCM Software manual and shared library in ~dist~ after compiling.
diff --git a/homeworks/hw-5.pdf b/homeworks/hw-5.pdf
new file mode 100644
index 0000000..ce43787
Binary files /dev/null and b/homeworks/hw-5.pdf differ
diff --git a/homeworks/hw-5.tex b/homeworks/hw-5.tex
new file mode 100644
index 0000000..8b9d24b
--- /dev/null
+++ b/homeworks/hw-5.tex
@@ -0,0 +1,95 @@
+% Created 2023-10-30 Mon 19:05
+% Intended LaTeX compiler: pdflatex
+\documentclass[11pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage{graphicx}
+\usepackage{longtable}
+\usepackage{wrapfig}
+\usepackage{rotating}
+\usepackage[normalem]{ulem}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{capt-of}
+\usepackage{hyperref}
+\notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+\author{Elizabeth Hunt}
+\date{\today}
+\title{Homework 5}
+\hypersetup{
+ pdfauthor={Elizabeth Hunt},
+ pdftitle={Homework 5},
+ pdfkeywords={},
+ pdfsubject={},
+ pdfcreator={Emacs 28.2 (Org mode 9.7-pre)}, 
+ pdflang={English}}
+\begin{document}
+
+\maketitle
+\setlength\parindent{0pt}
+
+\section{Question One}
+\label{sec:org88abf18}
+See LIZFCM \(\rightarrow\) Matrix Routines \(\rightarrow\) \texttt{lu decomp} \& \texttt{bsubst}.
+
+The test \texttt{UTEST(matrix, lu\_decomp)} is a unit test for the \texttt{lu\_decomp} routine,
+and \texttt{UTEST(matrix, bsubst)} verifies back substitution on an upper triangular
+3 \texttimes{} 3 matrix with a known solution that can be verified manually.
+
+Both can be found in \texttt{tests/matrix.t.c}.
+
+\section{Question Two}
+\label{sec:org098a7f1}
+Unless the following are met, the resulting solution will be garbage.
+
+\begin{enumerate}
+\item The matrix \(U\) must be not be singular.
+\item \(U\) must be square (or it will fail the \texttt{assert}).
+\item The system created by \(Ux = b\) must be consistent.
+\item \(U\) is in (obviously) upper-triangular form.
+\end{enumerate}
+
+Thus, the actual calculation performing the \(LU\) decomposition
+(in \texttt{lu\_decomp}) does a sanity
+check for 1-3 will fail an assert, should a point along the diagonal (pivot) be
+zero, or the matrix be non-factorable.
+
+\section{Question Three}
+\label{sec:org40d5983}
+See LIZFCM \(\rightarrow\) Matrix Routines \(\rightarrow\) \texttt{fsubst}.
+
+\texttt{UTEST(matrix, fsubst)} verifies forward substitution on a lower triangular 3 \texttimes{} 3
+matrix with a known solution that can be verified manually.
+
+\section{Question Four}
+\label{sec:orgf7d23bb}
+
+See LIZFCM \(\rightarrow\) Matrix Routines \(\rightarrow\) \texttt{gaussian\_elimination} and \texttt{solve\_gaussian\_elimination}.
+
+\section{Question Five}
+\label{sec:org54e966c}
+See LIZFCM \(\rightarrow\) Matrix Routines \(\rightarrow\) \texttt{m\_dot\_v}, and the \texttt{UTEST(matrix, m\_dot\_v)} in
+\texttt{tests/matrix.t.c}.
+
+\section{Question Six}
+\label{sec:org413b527}
+See \texttt{UTEST(matrix, solve\_gaussian\_elimination)} in \texttt{tests/matrix.t.c}, which generates a diagonally dominant 10 \texttimes{} 10 matrix
+and shows that the solution is consistent with the initial matrix, according to the steps given. Then,
+we do a dot product between each row of the diagonally dominant matrix and the solution vector to ensure
+it is near equivalent to the input vector.
+
+\section{Question Seven}
+\label{sec:orgd3d7443}
+See \texttt{UTEST(matrix, solve\_matrix\_lu\_bsubst)} which does the same test in Question Six with the solution according to
+\texttt{solve\_matrix\_lu\_bsubst} as shown in the Software Manual.
+
+\section{Question Eight}
+\label{sec:orgf8ac9bf}
+No, since the time complexity for Gaussian Elimination is always less than that of the LU factorization solution by \(O(n^2)\) operations
+(in LU factorization we perform both backwards and forwards substitutions proceeding the LU decomp, in Gaussian Elimination we only need
+back substitution).
+
+\section{Question Nine, Ten}
+\label{sec:orgb270171}
+See LIZFCM Software manual and shared library in \texttt{dist} after compiling.
+\end{document}
\ No newline at end of file
diff --git a/inc/lizfcm.h b/inc/lizfcm.h
index 12c1278..24b7fa9 100644
--- a/inc/lizfcm.h
+++ b/inc/lizfcm.h
@@ -25,6 +25,8 @@ extern double l2_distance(Array_double *v1, Array_double *v2);
 extern double l1_distance(Array_double *v1, Array_double *v2);
 extern double linf_distance(Array_double *v1, Array_double *v2);
 extern Array_double *copy_vector(Array_double *v1);
+extern Array_double *add_element(Array_double *v, double x);
+extern Array_double *slice_element(Array_double *v, size_t x);
 extern void free_vector(Array_double *v);
 extern void format_vector_into(Array_double *v, char *s);
 extern int vector_equal(Array_double *a, Array_double *b);
@@ -33,11 +35,15 @@ extern Matrix_double *put_identity_diagonal(Matrix_double *m);
 extern Matrix_double **lu_decomp(Matrix_double *m);
 extern Array_double *bsubst(Matrix_double *u, Array_double *b);
 extern Array_double *fsubst(Matrix_double *l, Array_double *b);
-extern Array_double *solve_matrix(Matrix_double *m, Array_double *b);
+extern Array_double *solve_matrix_lu_bsubst(Matrix_double *m, Array_double *b);
+extern Matrix_double *gaussian_elimination(Matrix_double *m);
+extern Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b);
 extern Array_double *m_dot_v(Matrix_double *m, Array_double *v);
 extern Matrix_double *m_dot_m(Matrix_double *a, Matrix_double *b);
 extern Array_double *col_v(Matrix_double *m, size_t x);
 extern Matrix_double *copy_matrix(Matrix_double *m);
+extern Matrix_double *add_column(Matrix_double *m, Array_double *col);
+extern Matrix_double *slice_column(Matrix_double *m, size_t col);
 extern void free_matrix(Matrix_double *m);
 extern void format_matrix_into(Matrix_double *m, char *s);
 extern int matrix_equal(Matrix_double *a, Matrix_double *b);
diff --git a/inc/macros.h b/inc/macros.h
index eab1b41..d081869 100644
--- a/inc/macros.h
+++ b/inc/macros.h
@@ -52,7 +52,7 @@
 #define c_max(x, y) (((x) >= (y)) ? (x) : (y))
 #define c_min(x, y) (((x) <= (y)) ? (x) : (y))
 
-#define true 1;
-#define false 0;
+#define true 1
+#define false 0
 
 #endif // MACROS_H
diff --git a/notes/Oct-27.org b/notes/Oct-27.org
new file mode 100644
index 0000000..6d23576
--- /dev/null
+++ b/notes/Oct-27.org
@@ -0,0 +1,26 @@
+Use a bisection criterion for a start
+
+Hybrid Method: combine Bisection and Higher Order Method:
+- Newton's Method
+- Secant Method (Newton's method with secant approx.)
+
+
+#+BEGIN_SRC c
+fa = f(a)
+fb = f(b)
+if (fa * fb >= 0) return
+
+error = 10 * tol
+iter = 0
+
+while (error > tol && iter < maxiter) {
+x0 = 0.5 * (a + b)
+x1 = x0 - f(x0) / f'(x0)
+if (abs(x1 - x0) > 0.5 * (b - a)) {
+// do bisection
+} else{
+// do newton's method
+}
+}
+#+END_SRC
+
diff --git a/notes/Oct-30.org b/notes/Oct-30.org
new file mode 100644
index 0000000..7d6ee03
--- /dev/null
+++ b/notes/Oct-30.org
@@ -0,0 +1,34 @@
+* Power Method for computing the largest eigenvalue of a square matrix
+
+An eigenvector, v \in R^n is a nonzero vector such that for some number, \lambda \in C, Av = \lambda v
+\Rightarrow || v || = 1
+
+
+Suppose we start with some vector v and assume, v = \alpha_0 v_0 + \alpha_1 v_1 + \cdots + \alpha_n v_n, where {v_1, \cdots, v_n}
+are the eigenvectors of A. Assume {v_1, \cdots, v_n} is a basis for R^n
+
+We can order the eigenvalues such that \lambda_1 \ge \lambda_2 \ge \lambda_3 \ge \cdots \ge \lambda_n
+
+Compute u = Av
+= A(\alpha_1 v_1 + \cdots + \alpha_n v_n)
+= \alpha_1 Av_1 + A(\cdots) + \alpha_n A v_n
+= \alpha_1 \lambda_1 v_1 + \alpha_2 \lambda_2 v_2 + \cdots + \alpha_n \lambda_n v_n
+
+w = A (Av)
+= \alpha_1 \lambda_1^2 v_1 + \alpha_2 \lambda_2^2 v_2 + \cdots + \alpha_n \lambda_n^2 v_n
+
+Thus,
+A^k v = \alpha_1 \lambda_1^k v_1 + \alpha_2 \lambda_2^k v_2 + \cdots + \alpha_n \lambda_n^k v_n
+= \lambda_1^k ( \alpha_1 v_1 + \alpha_2 \frac{\lambda_2^k}{\lambda_1^k} v_2 + \cdots + \alpha_n \frac{\lambda_3^k}{\lambda_1^k} v_n)
+
+As k \rightarrow \infty
+A^k v = \lambda_1^k (\alpha_1 v_1) + \text{negligble terms}
+
+Algorithm:
+v \ne 0 with v \in R^n
+y = Av = \alpha_1 v_1 + \cdots + \alpha_n v_n
+
+w = \frac{1}{||y||} \cdot y
+
+Rayleigh Quotient:
+If $v$ is an eigenvector of A with eigenvalue \lambda then \frac{v^T A v}{v^T v} = \lambda
diff --git a/src/matrix.c b/src/matrix.c
index 22dd171..0891734 100644
--- a/src/matrix.c
+++ b/src/matrix.c
@@ -1,5 +1,6 @@
 #include "lizfcm.h"
 #include <assert.h>
+#include <math.h>
 #include <stdio.h>
 #include <string.h>
 
@@ -71,7 +72,7 @@ Matrix_double **lu_decomp(Matrix_double *m) {
   for (size_t y = 0; y < m->rows; y++) {
     if (u->data[y]->data[y] == 0) {
       printf("ERROR: a pivot is zero in given matrix\n");
-      exit(-1);
+      assert(false);
     }
   }
 
@@ -82,7 +83,7 @@ Matrix_double **lu_decomp(Matrix_double *m) {
 
         if (denom == 0) {
           printf("ERROR: non-factorable matrix\n");
-          exit(-1);
+          assert(false);
         }
 
         double factor = -(u->data[y]->data[x] / denom);
@@ -129,7 +130,7 @@ Array_double *fsubst(Matrix_double *l, Array_double *b) {
   return x;
 }
 
-Array_double *solve_matrix(Matrix_double *m, Array_double *b) {
+Array_double *solve_matrix_lu_bsubst(Matrix_double *m, Array_double *b) {
   assert(b->size == m->rows);
   assert(m->rows == m->cols);
 
@@ -144,10 +145,99 @@ Array_double *solve_matrix(Matrix_double *m, Array_double *b) {
 
   free_matrix(u);
   free_matrix(l);
+  free(u_l);
 
   return x;
 }
 
+Matrix_double *gaussian_elimination(Matrix_double *m) {
+  uint64_t h = 0;
+  uint64_t k = 0;
+
+  Matrix_double *m_cp = copy_matrix(m);
+
+  while (h < m_cp->rows && k < m_cp->cols) {
+    uint64_t max_row = 0;
+    double total_max = 0.0;
+
+    for (uint64_t row = h; row < m_cp->rows; row++) {
+      double this_max = c_max(fabs(m_cp->data[row]->data[k]), total_max);
+      if (c_max(this_max, total_max) == this_max) {
+        max_row = row;
+      }
+    }
+
+    if (max_row == 0) {
+      k++;
+      continue;
+    }
+
+    Array_double *swp = m_cp->data[max_row];
+    m_cp->data[max_row] = m_cp->data[h];
+    m_cp->data[h] = swp;
+
+    for (uint64_t row = h + 1; row < m_cp->rows; row++) {
+      double factor = m_cp->data[row]->data[k] / m_cp->data[h]->data[k];
+      m_cp->data[row]->data[k] = 0.0;
+
+      for (uint64_t col = k + 1; col < m_cp->cols; col++) {
+        m_cp->data[row]->data[col] -= m_cp->data[h]->data[col] * factor;
+      }
+    }
+
+    h++;
+    k++;
+  }
+
+  return m_cp;
+}
+
+Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b) {
+  assert(b->size == m->rows);
+  assert(m->rows == m->cols);
+
+  Matrix_double *m_augment_b = add_column(m, b);
+  Matrix_double *eliminated = gaussian_elimination(m_augment_b);
+
+  Array_double *b_gauss = col_v(eliminated, m->cols);
+  Matrix_double *u = slice_column(eliminated, m->rows);
+
+  Array_double *solution = bsubst(u, b_gauss);
+
+  free_matrix(m_augment_b);
+  free_matrix(eliminated);
+  free_matrix(u);
+  free_vector(b_gauss);
+
+  return solution;
+}
+
+Matrix_double *slice_column(Matrix_double *m, size_t x) {
+  Matrix_double *sliced = copy_matrix(m);
+
+  for (size_t row = 0; row < m->rows; row++) {
+    Array_double *old_row = sliced->data[row];
+    sliced->data[row] = slice_element(old_row, x);
+    free_vector(old_row);
+  }
+  sliced->cols--;
+
+  return sliced;
+}
+
+Matrix_double *add_column(Matrix_double *m, Array_double *v) {
+  Matrix_double *pushed = copy_matrix(m);
+
+  for (size_t row = 0; row < m->rows; row++) {
+    Array_double *old_row = pushed->data[row];
+    pushed->data[row] = add_element(old_row, v->data[row]);
+    free_vector(old_row);
+  }
+
+  pushed->cols++;
+  return pushed;
+}
+
 void free_matrix(Matrix_double *m) {
   for (size_t y = 0; y < m->rows; ++y)
     free_vector(m->data[y]);
diff --git a/src/vector.c b/src/vector.c
index 3e4f62d..1b3e0b0 100644
--- a/src/vector.c
+++ b/src/vector.c
@@ -88,6 +88,21 @@ Array_double *copy_vector(Array_double *v) {
   return copy;
 }
 
+Array_double *add_element(Array_double *v, double x) {
+  Array_double *pushed = InitArrayWithSize(double, v->size + 1, 0.0);
+  for (size_t i = 0; i < v->size; ++i)
+    pushed->data[i] = v->data[i];
+  pushed->data[v->size] = x;
+  return pushed;
+}
+
+Array_double *slice_element(Array_double *v, size_t x) {
+  Array_double *sliced = InitArrayWithSize(double, v->size - 1, 0.0);
+  for (size_t i = 0; i < v->size - 1; ++i)
+    sliced->data[i] = i >= x ? v->data[i + 1] : v->data[i];
+  return sliced;
+}
+
 void free_vector(Array_double *v) {
   free(v->data);
   free(v);
diff --git a/test/matrix.t.c b/test/matrix.t.c
index 5386635..1c72b85 100644
--- a/test/matrix.t.c
+++ b/test/matrix.t.c
@@ -7,6 +7,38 @@ UTEST(matrix, free) {
   EXPECT_NE(data_addr, (uint64_t)(m->data));
 }
 
+UTEST(matrix, add_column) {
+  Matrix_double *m = InitMatrixWithSize(double, 5, 5, 0.0);
+  Array_double *col = InitArray(double, {1.0, 2.0, 3.0, 4.0, 5.0});
+  Matrix_double *new_m = add_column(m, col);
+
+  for (size_t row = 0; row < m->rows; row++)
+    EXPECT_EQ(new_m->data[row]->data[m->cols], col->data[row]);
+  EXPECT_EQ(new_m->cols, m->cols + 1);
+
+  free_matrix(m);
+  free_matrix(new_m);
+  free_vector(col);
+}
+
+UTEST(matrix, slice_column) {
+  size_t slice = 1;
+
+  Matrix_double *m = InitMatrixWithSize(double, 5, 5, 1.0 * (rand() % 10));
+  Matrix_double *new_m = slice_column(m, slice);
+
+  for (size_t row = 0; row < m->rows; row++) {
+    Array_double *sliced_row = slice_element(m->data[row], slice);
+
+    EXPECT_TRUE(vector_equal(new_m->data[row], sliced_row));
+    free_vector(sliced_row);
+  }
+  EXPECT_EQ(new_m->cols, m->cols - 1);
+
+  free_matrix(m);
+  free_matrix(new_m);
+}
+
 UTEST(matrix, put_identity_diagonal) {
   Matrix_double *m = InitMatrixWithSize(double, 8, 8, 0.0);
   Matrix_double *ident = put_identity_diagonal(m);
@@ -47,11 +79,53 @@ UTEST(matrix, m_dot_v) {
   free_vector(dotted);
 }
 
+UTEST(matrix, bsubst) {
+  Matrix_double *u = InitMatrixWithSize(double, 3, 3, 0.0);
+  u->data[0]->data[0] = 1.0;
+  u->data[0]->data[1] = 2.0;
+  u->data[0]->data[2] = 3.0;
+  u->data[1]->data[1] = 4.0;
+  u->data[1]->data[2] = 5.0;
+  u->data[2]->data[2] = 6.0;
+
+  Array_double *b = InitArray(double, {14.0, 29.0, 30.0});
+
+  Array_double *solution = bsubst(u, b);
+  EXPECT_NEAR(solution->data[0], -3.0, 0.0001);
+  EXPECT_NEAR(solution->data[1], 1.0, 0.0001);
+  EXPECT_NEAR(solution->data[2], 5.0, 0.0001);
+
+  free_matrix(u);
+  free_vector(b);
+  free_vector(solution);
+}
+
+UTEST(matrix, fsubst) {
+  Matrix_double *l = InitMatrixWithSize(double, 3, 3, 0.0);
+  l->data[0]->data[0] = 1.0;
+  l->data[1]->data[0] = 2.0;
+  l->data[1]->data[1] = 3.0;
+  l->data[2]->data[0] = 4.0;
+  l->data[2]->data[1] = 5.0;
+  l->data[2]->data[2] = 6.0;
+
+  Array_double *b = InitArray(double, {14.0, 13.0, 32.0});
+
+  Array_double *solution = fsubst(l, b);
+  EXPECT_NEAR(solution->data[0], 14.0, 0.0001);
+  EXPECT_NEAR(solution->data[1], -5.0, 0.0001);
+  EXPECT_NEAR(solution->data[2], 0.16667, 0.0001);
+
+  free_matrix(l);
+  free_vector(b);
+  free_vector(solution);
+}
+
 UTEST(matrix, lu_decomp) {
-  Matrix_double *m = InitMatrixWithSize(double, 8, 8, 0.0);
+  Matrix_double *m = InitMatrixWithSize(double, 10, 10, 0.0);
   for (size_t y = 0; y < m->rows; ++y) {
     for (size_t x = 0; x < m->cols; ++x)
-      m->data[y]->data[x] = x == y ? 5.0 : (5.0 - rand() % 10 + 1);
+      m->data[y]->data[x] = x == y ? 20.0 : (100.0 - rand() % 100) / 100.0;
   }
 
   Matrix_double **ul = lu_decomp(m);
@@ -75,15 +149,40 @@ UTEST(matrix, lu_decomp) {
   free(ul);
 }
 
-UTEST(matrix, solve_matrix) {
-  Matrix_double *m = InitMatrixWithSize(double, 8, 8, 0.0);
+UTEST(matrix, solve_gaussian_elimination) {
+  Matrix_double *m = InitMatrixWithSize(double, 10, 10, 0.0);
   for (size_t y = 0; y < m->rows; ++y) {
     for (size_t x = 0; x < m->cols; ++x)
-      m->data[y]->data[x] = x == y ? 10.0 : (5.0 - rand() % 10 + 1);
+      m->data[y]->data[x] = x == y ? 20.0 : (100.0 - rand() % 100) / 100.0;
   }
 
-  Array_double *b = InitArray(double, {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0});
-  Array_double *solution = solve_matrix(m, b);
+  Array_double *b_1 = InitArrayWithSize(double, m->rows, 1.0);
+  Array_double *b = m_dot_v(m, b_1);
+
+  Array_double *solution = solve_matrix_gaussian(m, b);
+
+  for (size_t y = 0; y < m->rows; y++) {
+    double dot = v_dot_v(m->data[y], solution);
+    EXPECT_NEAR(b->data[y], dot, 0.0001);
+  }
+
+  free_vector(b_1);
+  free_matrix(m);
+  free_vector(b);
+  free_vector(solution);
+}
+
+UTEST(matrix, solve_matrix_lu_bsubst) {
+  Matrix_double *m = InitMatrixWithSize(double, 10, 10, 0.0);
+  for (size_t y = 0; y < m->rows; ++y) {
+    for (size_t x = 0; x < m->cols; ++x)
+      m->data[y]->data[x] = x == y ? 20.0 : (100.0 - rand() % 100) / 100.0;
+  }
+
+  Array_double *b_1 = InitArrayWithSize(double, m->rows, 1.0);
+  Array_double *b = m_dot_v(m, b_1);
+
+  Array_double *solution = solve_matrix_lu_bsubst(m, b);
 
   for (size_t y = 0; y < m->rows; y++) {
     double dot = v_dot_v(m->data[y], solution);
@@ -92,6 +191,7 @@ UTEST(matrix, solve_matrix) {
 
   free_matrix(m);
   free_vector(b);
+  free_vector(b_1);
   free_vector(solution);
 }
 
diff --git a/test/vector.t.c b/test/vector.t.c
index 5dc8ba9..4811113 100644
--- a/test/vector.t.c
+++ b/test/vector.t.c
@@ -10,6 +10,28 @@ UTEST(vector, copy_vector) {
   free_vector(w);
 }
 
+UTEST(vector, add_element) {
+  Array_double *v = InitArray(double, {3, 1, -4});
+  Array_double *w = add_element(v, -2);
+  Array_double *w_expect = InitArray(double, {3, 1, -4, -2});
+  EXPECT_TRUE(vector_equal(w, w_expect));
+
+  free_vector(v);
+  free_vector(w);
+  free_vector(w_expect);
+}
+
+UTEST(vector, slice_element) {
+  Array_double *v = InitArray(double, {3, 1, -4});
+  Array_double *w = slice_element(v, 1);
+  Array_double *w_expect = InitArray(double, {3, -4});
+  EXPECT_TRUE(vector_equal(w, w_expect));
+
+  free_vector(v);
+  free_vector(w);
+  free_vector(w_expect);
+}
+
 UTEST(vector, free_vector) {
   Array_double *v = InitArray(double, {3, 1, -4});
   uint64_t arr_addr = (uint64_t)v->data;