diff --git a/code/4-matrix.R b/code/4-matrix.R
index 3cf53fe..65715e4 100644
--- a/code/4-matrix.R
+++ b/code/4-matrix.R
@@ -30,4 +30,6 @@ while (steps < 100) {
 
 cat(sprintf("Took %d steps\n", steps))
 cat(sprintf("Final matrix:\n"))
+step_matrix <- apply(step_matrix, 1, function(cell) round(cell, 4))
 print(step_matrix)
+write.table(step_matrix, "output/4-matrix.csv", sep = "\t", col.names = FALSE, row.names = FALSE)
diff --git a/code/4-sim.R b/code/4-sim.R
index 68bea82..adae644 100644
--- a/code/4-sim.R
+++ b/code/4-sim.R
@@ -25,8 +25,7 @@ for (round_idx in 1:rounds) {
   cur_state <- next_state
 }
 
-for (state in 1:3) {
-  state_count <- length(states[states == state])
-  state_prob <- state_count / rounds
-  cat(sprintf("State %d prob: %f\n", state, state_prob))
-}
+probs <- sapply(1:3, function(state) length(states[states == state]) / length(states))
+probs <- sapply(probs, function(x) round(x, 4))
+print(probs)
+write.table(probs, "output/4-sim.csv", sep = "\t", col.names = FALSE, row.names = FALSE)
diff --git a/code/4-solve.R b/code/4-solve.R
index 78227f5..c25482a 100644
--- a/code/4-solve.R
+++ b/code/4-solve.R
@@ -11,16 +11,7 @@ prob_matrix <- matrix(
 )
 
 # Build the coefficient matrix and constants
-coef_matrix <- matrix(
-  c(
-    c(-prob_matrix[1, 2] - prob_matrix[1, 3], prob_matrix[2, 1], prob_matrix[3, 1]),
-    c(prob_matrix[1, 2], -prob_matrix[2, 1] - prob_matrix[2, 3], prob_matrix[3, 2]),
-    c(prob_matrix[1, 3], prob_matrix[2, 3], -prob_matrix[3, 1] - prob_matrix[3, 2])
-  ),
-  nrow = 3,
-  ncol = 3,
-  byrow = TRUE
-)
+coef_matrix <- t(prob_matrix) - diag(rowSums(prob_matrix))
 consts <- c(0, 0, 0)
 
 # Substitute one of the rows by the final `pi0 + pi1 + pi2 = 1` equation.
@@ -29,4 +20,6 @@ coef_matrix[subst_idx, ] <- c(1, 1, 1)
 consts[subst_idx] <- 1
 
 solution <- solve(coef_matrix, consts)
+solution <- sapply(solution, function(x) round(x, 4))
 print(solution)
+write.table(solution, "output/4-solve.csv", sep = "\t", col.names = FALSE, row.names = FALSE)
diff --git a/typst/exercises/2.typ b/typst/exercises/2.typ
index 3cca4ab..23de6e0 100644
--- a/typst/exercises/2.typ
+++ b/typst/exercises/2.typ
@@ -22,7 +22,7 @@
 #indent_par[Afterwards, we can divide each row by the number of occurrences in that row to obtain the transition probability matrix. The following code 1 is the code we developed to accomplish this:]
 
 #code_figure(
-  text(size: 0.88em, raw(read("/code/2.R"), lang: "R", block: true)),
+  text(size: 0.8em, raw(read("/code/2.R"), lang: "R", block: true)),
   caption: "Developed code",
 )
 
@@ -34,7 +34,6 @@
   columns: (1fr, 1fr),
   figure(
     pad(1em, text(size: 1.8em, math.mat(
-      delim: "[",
       gap: 1em,
       ..occur_matrix
     ))),
@@ -43,7 +42,6 @@
   ),
   figure(
     pad(1em, text(size: 1.8em, math.mat(
-      delim: "[",
       gap: 1em,
       ..prob_matrix
     ))),
diff --git a/typst/exercises/3.typ b/typst/exercises/3.typ
index a71415f..7caff31 100644
--- a/typst/exercises/3.typ
+++ b/typst/exercises/3.typ
@@ -41,3 +41,5 @@
 )
 
 ==== c.
+
+#pagebreak()
diff --git a/typst/exercises/4.typ b/typst/exercises/4.typ
index 2139c26..2c82ac8 100644
--- a/typst/exercises/4.typ
+++ b/typst/exercises/4.typ
@@ -1,5 +1,89 @@
-#import "/typst/util.typ" as util: indent_par
+#import "/typst/util.typ" as util: indent_par, code_figure
 
 ==== (i) Solving the balance equations
 
-#indent_par[]
+#indent_par[The following equations 1, 2, 3 and 4 are our balance equations:]
+
+$ pi_0 p_00 + pi_1 p_10 + pi_2 p_20 = pi_0 p_00 + pi_0 p_01 + pi_0 p_02 $
+$ pi_0 p_01 + pi_1 p_11 + pi_2 p_21 = pi_1 p_10 + pi_1 p_11 + pi_1 p_12 $
+$ pi_0 p_02 + pi_1 p_12 + pi_2 p_22 = pi_2 p_20 + pi_2 p_21 + pi_2 p_22 $
+$ pi_0 + pi_1 + pi_2 = 1 $
+
+#indent_par[Calling $P$ the matrix $mat(p_00, p_01, p_02; p_10, p_11, p_12; p_20, p_21, p_22)$, the 3 first balance equations can be expressed as the following equation 5:]
+
+$ P^T dot mat(sum_(j) p_0j, 0, 0; 0, sum_(j) p_1j, 0; 0, 0, sum_(j) p_2j;) = 0 $
+
+#indent_par[We can then model this in R using the following code 2:]
+
+#code_figure(
+  text(size: 0.8em, raw(read("/code/4-solve.R"), lang: "R", block: true)),
+  caption: "Code for solving the balance equations",
+)
+
+#indent_par[To solve the system, the 3 first balance equations aren't enough, we must substitute one of them for the 4th balance equations to ensure the system is solvable. In our case, we substituted the 3rd equation with the 4th. The following table 3 contains our results:]
+
+#let solution = csv("/output/4-solve.csv", delimiter: "\t")
+#figure(
+	pad(1em, text(size: 1.5em, math.mat(
+		gap: 1em,
+		..solution
+	))),
+	kind: table,
+	caption: "Limiting state probabilities"
+)
+
+==== (ii). Matrix multiplication
+
+#indent_par[The following code 3 contains our approach to obtain the limiting state probabilities via matrix multiplication.]
+
+#code_figure(
+  text(size: 0.8em, raw(read("/code/4-matrix.R"), lang: "R", block: true)),
+  caption: "Code using matrix multiplication",
+)
+
+#indent_par[We can simply raise the probability transition matrix by itself numerous times to get a matrix that has our limiting state probabilities in each row.]
+
+#indent_par[To decide how many times we need to raise the matrix, we keep raising it until the difference in each cell from the previous step is less than a certain $ε$. In our case, we used $ε = 10^(-9)$. This process took, in total 22 steps to achieve our desired $ε$, obtaining the limiting state probabilities in table 4:]
+
+#let solution = csv("/output/4-matrix.csv", delimiter: "\t")
+#figure(
+	pad(1em, text(size: 1.5em, math.mat(
+		gap: 1em,
+		..solution
+	))),
+	kind: table,
+	caption: "Limiting state probabilities"
+)
+
+#indent_par[The results we obtained are equal to the theoretical results calculated in the previous exercise.]
+
+#pagebreak()
+
+==== (iii). Simulation
+
+#indent_par[The following code 4 contains our approach to obtain the limiting state probabilities via simulation.]
+
+#code_figure(
+  text(size: 0.8em, raw(read("/code/4-sim.R"), lang: "R", block: true)),
+  caption: "Code using simulation",
+)
+
+#indent_par[We initialize our current state to 1, then for 100000 rounds, save the current state, calculate the next state and save it to the current.]
+
+#indent_par[In order to calculate the next state, we generate a uniformly random number in the $[0.0, 1.0]$ interval, and then choose the first index of the cumulative sum of the probabilities that is higher than the number we generated.]
+
+#indent_par[This works because by calculating the cumulative sum of the probabilities, we're calculating it's cumulative density function. Then by definition, finding the input for which this function has value $<= x$, for $x ∈ [0.0, 1.0]$ is equal to sampling the original distribution.]
+
+#indent_par[Finally, we get the limiting state probabilities by checking how many states there are out of all the ones we've visited for each state. The output of this process can be seen in table 5:]
+
+#let solution = csv("/output/4-sim.csv", delimiter: "\t")
+#figure(
+	pad(1em, text(size: 1.5em, math.mat(
+		gap: 1em,
+		..solution
+	))),
+	kind: table,
+	caption: "Limiting state probabilities"
+)
+
+#indent_par[The results aren't exactly equal to the theoretical results, but this is expected, given it's a simulation. However, they are very close to them.]
diff --git a/typst/main.typ b/typst/main.typ
index b2eaa55..03f307e 100644
--- a/typst/main.typ
+++ b/typst/main.typ
@@ -20,6 +20,13 @@
 #set par(
   justify: true,
 )
+#set math.equation(
+  numbering: "(1)",
+)
+#set math.mat(
+  delim: "[",
+
+)
 #show link: underline
 
 #include "cover.typ"
diff --git a/zbuild.yaml b/zbuild.yaml
index b2304f0..287c210 100644
--- a/zbuild.yaml
+++ b/zbuild.yaml
@@ -1,5 +1,13 @@
 ---
-default: [output/1.svg, rule: ex1, rule: ex2, rule: ex3_a, rule: ex3_b]
+default:
+  - output/1.svg
+  - rule: ex1
+  - rule: ex2
+  - rule: ex3_a
+  - rule: ex3_b
+  - rule: ex4_solve
+  - rule: ex4_matrix
+  - rule: ex4_sim
 
 rules:
   # Typst
@@ -75,3 +83,33 @@ rules:
     exec:
       - - Rscript
         - code/3b.R
+
+  # Exercise 4 solve
+  ex4_solve:
+    out:
+      - output/4-solve.csv
+    deps:
+      - code/4-solve.R
+    exec:
+      - - Rscript
+        - code/4-solve.R
+
+  # Exercise 4 matrix
+  ex4_matrix:
+    out:
+      - output/4-matrix.csv
+    deps:
+      - code/4-matrix.R
+    exec:
+      - - Rscript
+        - code/4-matrix.R
+
+  # Exercise 4 sim
+  ex4_sim:
+    out:
+      - output/4-sim.csv
+    deps:
+      - code/4-sim.R
+    exec:
+      - - Rscript
+        - code/4-sim.R