diff --git a/week3/NMM Random Systems_230302_141008.pdf b/week3/NMM Random Systems_230302_141008.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..7fcdf86785af4d954f996c6439e728fdf5f13f82
Binary files /dev/null and b/week3/NMM Random Systems_230302_141008.pdf differ
diff --git a/week3/README.MD b/week3/README.MD
new file mode 100644
index 0000000000000000000000000000000000000000..353d3867626168ea1b5c87f3bb855cbc39d9a296
--- /dev/null
+++ b/week3/README.MD
@@ -0,0 +1,7 @@
+# Week 3 - Random Systems
+
+[Here's the written portion](NMM_ODEs_230223_121852.pdf)
+
+Code is [here](random_walker.py).
+
+![](img/random_walker.jpg)
diff --git a/week3/random_walker.jpg b/week3/random_walker.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..7648dfbf0634e31a876e483b25e8c718e734f6c2
Binary files /dev/null and b/week3/random_walker.jpg differ
diff --git a/week3/random_walker.py b/week3/random_walker.py
new file mode 100644
index 0000000000000000000000000000000000000000..24e797b4c15a206497a1032adb3a7ce239da5828
--- /dev/null
+++ b/week3/random_walker.py
@@ -0,0 +1,60 @@
+'''
+Part D:
+Write a program (including the random number generator) to plot the position as a function of time of a random walker
+in 1D that at each time step has an equal probability of making a step of ±1.
+Plot an ensemble of 10 trajectories, each 1000 points long, and overlay error bars of width 3sigma(t) on the plot.
+'''
+import numpy as np
+import time
+
+def build_LFSR(seed, bits_needed=1000):
+    depth  = 24
+    taps   = [1,2,7,24]
+    state = seed
+    output = []
+    for _ in range(bits_needed):
+        newbit = (state ^ (state >> taps[1]) ^ (state >> taps[2]) ^ (state >> taps[3])) & 1
+        state = (state >> 1) | (newbit << (depth - 1)) # shift everything over once and or our new bit in
+        output.append(newbit)
+    return output
+
+# Initialize a numpy array of zeros 10x1000
+walker_states = np.zeros((10, 1000))
+
+# create our seeds by grabbing a 16 bit number from the LSB of the current time in us
+# cheating and just adding a fixed delay to get different seeds
+seeds = []
+for i in range(10):
+    delay = .1
+    seeds.append(int(str(time.time())[-3:]))
+    time.sleep(delay)
+
+for i in range(10):
+    random_path = np.array(build_LFSR(seeds[i]))
+    random_path[(random_path == 0)] = -1 # turn our os to -1s so we can actually walk
+    for j in range(1000):
+        walker_states[i, j] = walker_states[i, j-1] + random_path[j]
+
+# Brownian motion is normally distributed with variance = t, so standard deviation will be sqrt(t), so let's get all of our timesteps now
+timestep = np.array(range(1000))
+
+import matplotlib.pyplot as plt
+plt.plot(walker_states.T) # transpose it so we can plot it
+
+# for the 3 sigma lines, we will plot the sqrt of the timestep and divide the result by 1.5 to plot above and below the mean (0)
+plt.plot(timestep, 1.5 * np.sqrt(timestep), color='black', alpha = 0.25)
+plt.plot(timestep, -1.5 * np.sqrt(timestep), color='black', alpha = 0.25)
+plt.show()
+
+'''
+Part E: What fraction of the trajectories should be contained in the error bars?
+
+# We can refer to a standard normal distribution table (alternative to integrating)
+Z-tables have the first two digits of Z-score in Y and the following 2 in X, so for 1.5 we look at the first entry of the +1.5 column
+and the first entry of the +1.5 row, which is 0.93319. (from wikipedia)
+https://en.wikipedia.org/wiki/Standard_normal_table
+
+Our actual answer is:
+'''
+probability_within = 0.93319 - (1 - 0.93319)
+print(f"{probability_within*100:.2f} % to be within 3 sigma at any given time.")