diff --git a/week5/img/laplace_2D_edges.gif b/week5/img/laplace_2D_edges.gif
index 9cabe0b91b8d20fa25c060277d01c5d06d6734f7..92f86d158e98b730a388a4ff7a0b9d7a8f9842a7 100644
Binary files a/week5/img/laplace_2D_edges.gif and b/week5/img/laplace_2D_edges.gif differ
diff --git a/week5/week5_4.py b/week5/week5_4.py
index d6b5505218aff88adc13225d68bf3ba5f5f85d0f..227abd8476bc54843437f26f3f70b18f79d5c9c8 100644
--- a/week5/week5_4.py
+++ b/week5/week5_4.py
@@ -57,8 +57,8 @@ def update(i):
     ax.plot_surface(X, Y, u, cmap = 'viridis', edgecolor = 'none')
 
 anim = animation.FuncAnimation(fig, update, total_time, interval=1000/fps, blit = False)
-# plt.show()
+plt.show()
 
-f = r"c://Users/david/Desktop/NMM/week5/img/laplace_2D_edges.gif" 
-writervideo = animation.FFMpegWriter(fps=fps/2, extra_args=['-vcodec', 'libx264']) 
-anim.save(f, writer='imagemagick', fps=fps/3, dpi=96)
+# f = r"c://Users/david/Desktop/NMM/week5/img/laplace_2D_edges.gif" 
+# writervideo = animation.FFMpegWriter(fps=fps/2, extra_args=['-vcodec', 'libx264']) 
+# anim.save(f, writer='imagemagick', fps=fps/3, dpi=96)
diff --git a/week6/README.MD b/week6/README.MD
new file mode 100644
index 0000000000000000000000000000000000000000..0e1b8b91f03a0ab9589d63988bb4bc028442c6e1
--- /dev/null
+++ b/week6/README.MD
@@ -0,0 +1,3 @@
+# Week 6 - Finite Elements
+
+## 10.2
\ No newline at end of file
diff --git a/week6/week6_2.py b/week6/week6_2.py
new file mode 100644
index 0000000000000000000000000000000000000000..a049a74dffad1a6ba959594ed0e8575f16a66dce
--- /dev/null
+++ b/week6/week6_2.py
@@ -0,0 +1,121 @@
+import math
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+from matplotlib.animation import FuncAnimation, writers
+
+'''
+Adapted from: https://polymerfem.com/full-finite-element-solver-in-100-lines-of-python/
+6.2)
+Model the bending of a beam (equation 10.29) under an applied load. Use Hermite polynomial interpolation,
+and boundary conditions fixing the displacement and slope at one end, and applying a force at the other end.
+V = int(0, L)[0.5*EI * (d^2u/dx^2)^2 - u(x) * f(x)]dx
+'''
+
+# Supporting functions
+def shape(xi):
+	x,y = tuple(xi)
+	N = [(1.0-x)*(1.0-y), (1.0+x)*(1.0-y), (1.0+x)*(1.0+y), (1.0-x)*(1.0+y)]
+	return 0.25*np.array(N)
+
+def gradshape(xi):
+	x,y = tuple(xi)
+	dN = [[-(1.0-y),  (1.0-y), (1.0+y), -(1.0+y)],
+		  [-(1.0-x), -(1.0+x), (1.0+x),  (1.0-x)]]
+	return 0.25*np.array(dN)
+
+# Input
+mesh_ex = 49
+mesh_ey = 9
+mesh_lx = 50.0
+mesh_ly = 10.0
+# Derived
+mesh_nx      = mesh_ex + 1
+mesh_ny      = mesh_ey + 1
+num_nodes    = mesh_nx * mesh_ny
+num_elements = mesh_ex * mesh_ey
+mesh_hx      = mesh_lx / mesh_ex
+mesh_hy      = mesh_ly / mesh_ey
+nodes = [] # make our mesh nodes
+for y in np.linspace(0.0, mesh_ly, mesh_ny):
+	for x in np.linspace(0.0, mesh_lx, mesh_nx):
+		nodes.append([x,y]) # nodes is an array of X and Y coordinates
+nodes = np.array(nodes)
+
+conn = []
+for j in range(mesh_ey):
+	for i in range(mesh_ex):
+		n0 = i + j*mesh_nx
+		conn.append([n0, n0 + 1, n0 + 1 + mesh_nx, n0 + mesh_nx])
+
+# Material model - Plane strain
+E = 100.0
+v = 0.48
+C = E/(1.0+v)/(1.0-2.0*v) * np.array([[1.0-v,     v,     0.0],
+								      [    v, 1.0-v,     0.0],
+								      [  0.0,   0.0,   0.5-v]])
+
+# Create global stiffness matrix
+K = np.zeros((2*num_nodes, 2*num_nodes))
+q4 = [[x/math.sqrt(3.0),y/math.sqrt(3.0)] for y in [-1.0,1.0] for x in [-1.0,1.0]]
+B = np.zeros((3,8))
+for c in conn:
+	xIe = nodes[c,:]
+	Ke = np.zeros((8,8))
+	for q in q4:
+		dN = gradshape(q)
+		J  = np.dot(dN, xIe).T
+		dN = np.dot(np.linalg.inv(J), dN)
+		B[0,0::2] = dN[0,:]
+		B[1,1::2] = dN[1,:]
+		B[2,0::2] = dN[1,:]
+		B[2,1::2] = dN[0,:]
+		Ke += np.dot(np.dot(B.T,C),B) * np.linalg.det(J)
+	for i,I in enumerate(c):
+		for j,J in enumerate(c):
+			K[2*I, 2*J]     += Ke[2*i, 2*j]
+			K[2*I+1, 2*J]   += Ke[2*i+1, 2*j]
+			K[2*I+1, 2*J+1] += Ke[2*i+1, 2*j+1]
+			K[2*I, 2*J+1]   += Ke[2*i, 2*j+1]
+
+# Assign nodal forces and boundary conditions
+f = np.zeros((2*num_nodes))
+for i in range(num_nodes):
+	if nodes[i,0] == 0.0: # if we are at the far left, fix the displacement?
+		K[2*i, :]       = 0.0
+		K[2*i+1, :]     = 0.0
+		K[2*i, 2*i]     = 1.0
+		K[2*i+1, 2*i+1] = 1.0
+	if nodes[i,0] == mesh_lx: # if we are at the far right
+		y = nodes[i,1]      # get the y coordinate
+		f[2*i+1] = 0.5    # apply a force
+		if y == 0.0 or y == mesh_ly: # at the edges make the force half
+			f[2*i+1] *= 0.5
+
+# Solving linear system
+u = np.linalg.solve(K, f)
+print('max u=', max(u))
+
+# Plotting
+ux = np.reshape(u[0::2], (mesh_ny,mesh_nx))
+uy = np.reshape(u[1::2], (mesh_ny,mesh_nx))
+xvec = []
+yvec = []
+res  = []
+for i in range(mesh_nx):
+    for j in range(mesh_ny):
+        xvec.append(i*mesh_hx + ux[j,i])
+        yvec.append(j*mesh_hy + uy[j,i])
+        res.append(uy[j,i])
+# t = plt.tricontourf(xvec, yvec, res, levels=14, cmap=plt.cm.jet)
+plt.scatter(xvec, yvec, marker='o', c='black', s=2)
+# plt.grid()
+# plt.colorbar(t)
+plt.axis('equal')
+plt.show()
+
+# my_dpi = 96
+# f = r"c://Users/david/Desktop/NMM/week5/img/wave_no_damp.gif" 
+# fps2 = 20
+# writervideo = animation.FFMpegWriter(fps=fps2, extra_args=['-vcodec', 'libx264']) 
+# anim.save(f, writer='imagemagick', fps=fps2, dpi=my_dpi)
\ No newline at end of file
diff --git a/week7/README.MD b/week7/README.MD
new file mode 100644
index 0000000000000000000000000000000000000000..ec7a6e551f132c8b2faa91e4012a9021420a7b25
--- /dev/null
+++ b/week7/README.MD
@@ -0,0 +1,16 @@
+# Week 7 - Discrete Elements
+
+## 9_1
+
+<img src="img/not_a_beam.gif" height="600">
+<img src="img/wrong_fixed_beam.gif" height="300">
+
+![](img/beam_bending_large.mp4)
+
+## 9_2
+
+![](img/falling_objects.mp4)
+
+![](img/solid_fall_instability.mp4)
+
+[Talk by a Noita developer on their physics engine](https://www.youtube.com/watch?v=prXuyMCgbTc)
\ No newline at end of file
diff --git a/week7/figures/as.mp4 b/week7/figures/as.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..9d3959ea780c6d2d3b37fe970933a18a4297f2e1
Binary files /dev/null and b/week7/figures/as.mp4 differ
diff --git a/week7/figures/beam_bending.mp4 b/week7/figures/beam_bending.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..47f3ff51ce4a4c1162c0049498e9cd6803eeedbd
Binary files /dev/null and b/week7/figures/beam_bending.mp4 differ
diff --git a/week7/figures/beam_bending_large.mp4 b/week7/figures/beam_bending_large.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..aad3a816fb1c5187404a08afcfb65cab934a3770
Binary files /dev/null and b/week7/figures/beam_bending_large.mp4 differ
diff --git a/week7/figures/brittle_fall.mp4 b/week7/figures/brittle_fall.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..8cb9d7d7a7fa818f461557d8784c21a9630c8478
Binary files /dev/null and b/week7/figures/brittle_fall.mp4 differ
diff --git a/week7/figures/deformable_fall.mp4 b/week7/figures/deformable_fall.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..e7508b930ed3c3332ff325f7b7ab77a3f2c31812
Binary files /dev/null and b/week7/figures/deformable_fall.mp4 differ
diff --git a/week7/figures/falling_objects.mp4 b/week7/figures/falling_objects.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..73b145aa6c920407d9ac3ff02e8c19418b9cc125
Binary files /dev/null and b/week7/figures/falling_objects.mp4 differ
diff --git a/week7/figures/liquid_fall.mp4 b/week7/figures/liquid_fall.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..739ec827516407e4738b4cb49539e9918ce0752c
Binary files /dev/null and b/week7/figures/liquid_fall.mp4 differ
diff --git a/week7/figures/my_people.mp4 b/week7/figures/my_people.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..d1fa8533318749d258992a3c9ba54cc26e073737
Binary files /dev/null and b/week7/figures/my_people.mp4 differ
diff --git a/week7/figures/not_a_beam.gif b/week7/figures/not_a_beam.gif
new file mode 100644
index 0000000000000000000000000000000000000000..4c17d1ca522d072df0e82b2157672922ce258685
Binary files /dev/null and b/week7/figures/not_a_beam.gif differ
diff --git a/week7/figures/rigid_fall.mp4 b/week7/figures/rigid_fall.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..ffe600ec9d57dd0c6505e639f3517bed70769d87
Binary files /dev/null and b/week7/figures/rigid_fall.mp4 differ
diff --git a/week7/figures/solid_fall.mp4 b/week7/figures/solid_fall.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..2c0760a86a5aca733b4817a8c66cf1868d35a5de
Binary files /dev/null and b/week7/figures/solid_fall.mp4 differ
diff --git a/week7/figures/solid_fall_instability.mp4 b/week7/figures/solid_fall_instability.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..2d0cb4eb0a3181caed61b8f8aeaa1818417dc682
Binary files /dev/null and b/week7/figures/solid_fall_instability.mp4 differ
diff --git a/week7/figures/viscous_fall.mp4 b/week7/figures/viscous_fall.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..c7419ef6faa7ad9ee4da99e68b92cd1bc3f8a070
Binary files /dev/null and b/week7/figures/viscous_fall.mp4 differ
diff --git a/week7/figures/wrong_fixed_beam.gif b/week7/figures/wrong_fixed_beam.gif
new file mode 100644
index 0000000000000000000000000000000000000000..89ed8b13db88381c199c762fba62ebf118ea1fab
Binary files /dev/null and b/week7/figures/wrong_fixed_beam.gif differ
diff --git a/week7/week7_9_1.py b/week7/week7_9_1.py
new file mode 100644
index 0000000000000000000000000000000000000000..63da7757660c3e90bc2d91b0b18a6a76aea19f30
--- /dev/null
+++ b/week7/week7_9_1.py
@@ -0,0 +1,128 @@
+import scipy.spatial as spatial
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+
+def calc_force(P1,P2):
+    # calculate the force between points P1 and P2 (both 2D numpy arrays)
+    # the neutral distance is the variable neutral_distance
+    dist = np.linalg.norm(P1 - P2) # we should have this from scipy already but recalculating for now...
+    return spring_constant * (dist - neutral_distance) * (P1 - P2) / dist
+
+def solve_forces(particles, velocities):
+    particles_tree = spatial.cKDTree(particles)
+    # For each point, calculate the force due to the nearest neighbors
+    F = np.zeros((x_num*y_num, 2))
+    for particle in range(x_num*y_num):
+        # return the index of all particles within radius of radius_check
+        neighbors = particles_tree.query_ball_point(particles[particle, :], radius_check)
+        # calculate the force on the particle due to the nearest neighbor
+        for neighbor in neighbors:
+            if neighbor != particle:
+                F[particle, :] -= calc_force(particles[particle, :], particles[neighbor, :])
+    F_damping = damping_constant * velocities  # damping forces
+    # add damping forces to the total force, ensuring it can't reverse the direction of the force
+    F[:, 0] -= np.sign(velocities[:, 0]) * np.minimum(np.abs(F_damping[:, 0]), np.abs(F[:, 0]))
+    F[:, 1] -= np.sign(velocities[:, 1]) * np.minimum(np.abs(F_damping[:, 1]), np.abs(F[:, 1]))
+    F[:, 1] -= gravity  # gravity
+    return F
+
+def solve_velocities(F, V_prev):
+    # update the velocity of each particle from V_prev according to the force acting on it
+    V = V_prev + dt * F / particle_mass
+    return V
+
+def solve_positions(particles, V):
+    # update the position of each particle from particles according to the velocity acting on it
+    particles_new = particles + dt * V
+    # CANTILEVERED BEAM - set all particles in the first column to have a position of zero
+    for particle_num in range(len(particles)):
+        # if the particle is in the first row:
+        if particle_num <= y_num:
+            # keep that particle in place
+            particles_new[particle_num, 0] = particles[particle_num, 0]
+            particles_new[particle_num, 1] = particles[particle_num, 1]
+    return particles_new
+
+def best_fit(particles, order):
+    a, b = np.polyfit(particles[:, 0], particles[:, 1], deg = order)
+    return np.array([a, b])
+
+# Constants:
+x_num = 100
+y_num = 30
+t_final = 1.5
+dt = 0.01
+mesh_spacing = 1
+neutral_distance = mesh_spacing
+radius_check = mesh_spacing*1.5
+spring_constant = 1000
+damping_constant = 10
+gravity = 0.05
+particle_mass = 1
+height_offset = 8
+
+sub_increments = 10
+
+# Create our triangular mesh of points
+particles = np.zeros((x_num*y_num, 2))
+velocities = np.zeros((x_num*y_num, 2))
+forces = np.zeros((x_num*y_num, 2))
+
+# create a mesh of particles, indexes start at point 0,0 and grow up the y-axis, then move to the next column and grow up that one
+for i in range(x_num):
+    for j in range(y_num):
+        particles[i*y_num+j, 0] = i
+        particles[i*y_num+j, 1] = j * np.sqrt(3) / 2 # height of equilateral triangle
+        if j % 2 == 1:  # add an offset to every other point
+            particles[i*y_num+j, 0] += 0.5
+
+# particles[-1, 0] -= mesh_spacing*.5 # offset a particle for chaos
+particles *= mesh_spacing # scale the mesh to the desired spacing
+particles[:,1] += mesh_spacing*height_offset
+
+# Create a scatterplot of X
+my_dpi = 96
+fig = plt.figure(figsize=(800/my_dpi, 400/my_dpi), dpi=my_dpi)
+ax = fig.add_axes([0, 0, 1, 1], frameon=False)
+# numpy array of colors for each particle in shades from red to blue
+colors = np.zeros((x_num*y_num, 4))
+colors[:, 0] = np.linspace(1, 0, x_num*y_num)
+colors[:, 2] = np.linspace(0, 1, x_num*y_num)
+colors[:, 3] = 1
+plt.axis('equal')   # ensure axes are equally scaled
+plt.plot([-mesh_spacing*3, x_num*mesh_spacing+mesh_spacing*3], [0, 0], color = 'k') # Draw a line representing the ground at y = 0;
+best_fit_0 = best_fit(particles, 1)
+static_line = ax.plot([0, x_num*mesh_spacing], [best_fit_0[0]*0 + best_fit_0[1], best_fit_0[0]*x_num*mesh_spacing + best_fit_0[1]], color='grey', linestyle='--', linewidth=1)
+fit_line = ax.plot([0, x_num*mesh_spacing], [best_fit_0[0]*0 + best_fit_0[1], best_fit_0[0]*x_num*mesh_spacing + best_fit_0[1]], color='black', linestyle='--', linewidth=1)
+plot = ax.scatter(particles[:, 0], particles[:, 1],
+                  marker = 'o',
+                  s = 20,
+                  color = colors,
+)
+# Optionally plot the intial conditions
+# plt.scatter(particles[:, 0], particles[:, 1], marker = 'o', s = 10, color = 'b')
+
+def animate(w):
+    global particles, velocities, forces
+    for i in range(sub_increments):
+        forces = solve_forces(particles, velocities)
+        velocities = solve_velocities(forces, velocities)
+        particles_new = solve_positions(particles, velocities)
+        x = particles_new[:, 0]
+        y = particles_new[:, 1]
+        particles = particles_new  
+    plot.set_offsets(np.c_[x,y])
+    best_fit_new = best_fit(particles, 1)
+    # update fit_line with the new best fit data
+    fit_line[0].set_data([0, x_num*mesh_spacing], [best_fit_new[0]*0 + best_fit_new[1], best_fit_new[0]*x_num*mesh_spacing + best_fit_new[1]])
+    # plt.plot([0, x_num*mesh_spacing], , color='steelblue', linestyle='--', linewidth=1)
+
+anim = animation.FuncAnimation(fig, animate, frames=int(t_final/dt))
+# plt.show()
+
+import matplotlib as mpl 
+mpl.rcParams['animation.ffmpeg_path'] = r'C:\ffmpeg_2023_03_30\bin\ffmpeg.exe'
+path = r"c://Users/david/Desktop/NMM/week7/figures/animation.mp4" 
+writervideo = animation.FFMpegWriter(fps=30, extra_args=['-vcodec', 'libx264']) 
+anim.save(path, writer=writervideo, dpi=my_dpi)
\ No newline at end of file
diff --git a/week7/week7_9_2.py b/week7/week7_9_2.py
new file mode 100644
index 0000000000000000000000000000000000000000..e578ced9478d515b653d9cde42de1322272a9284
--- /dev/null
+++ b/week7/week7_9_2.py
@@ -0,0 +1,159 @@
+import scipy.spatial as spatial
+from scipy import integrate as integrate
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+import time
+
+# import jit from numba
+# from numba import jit
+
+# @jit(nopython=True)
+def calc_force(P1, P2):
+    # calculate the force between points P1 and P2 (both 2D numpy arrays)
+    # the neutral distance is the variable neutral_distance
+    dist = np.linalg.norm(P1 - P2) # we should have this from scipy already but recalculating for now...
+    force = spring_constant * (dist - neutral_distance) * (P1 - P2) / dist
+    # repulsive forces only
+    # force = np.where(dist < neutral_distance, force, force*0.1)
+    return force
+
+# @jit(nopython=True)
+def solve_forces(particles, velocities):
+    particles_tree = spatial.cKDTree(particles)
+    # For each point, calculate the force due to the nearest neighbors
+    F = np.zeros((x_num*y_num, 2))
+    for particle in range(x_num*y_num):
+        # return the index of all particles within radius of radius_check
+        neighbors = particles_tree.query_ball_point(particles[particle, :], radius_check)
+        # Moving local damping here to use average velocity to avoid overapplying it
+        avg_vel = np.mean(velocities[neighbors, :], axis=0) # axis = 0 should average the columns
+        F_damping = damping_constant * (velocities[particle, :] - avg_vel)
+        # calculate the force on the particle due its neighbors
+        for neighbor in neighbors:
+            if neighbor != particle:
+                F[particle, :] -= calc_force(particles[particle, :], particles[neighbor, :])
+        F[particle, :] -= F_damping # leaving this out here for now
+    # Add some global damping
+    F -= damping_constant_global * velocities
+    F[:, 1] -= gravity  # gravity
+    return F
+
+# @jit(nopython=True)
+def solve_velocities(F, particles, dt):
+    # update the velocity of each particle from V_prev according to the force acting on it
+    V = dt * F / particle_mass
+    # if we contact the ground, reverse the y velocity and dampen it
+    V[:, 1] = np.where((particles[:, 1] < 0), -V[:, 1]*coef_of_restitution, V[:, 1]) # if we cross into the floor, reverse the y velocity and apply a damping factor
+    # apply friction to the ground
+    V[:, 0] = np.where((particles[:, 1] < 0), V[:, 0]*0.9, V[:, 0])
+    return V
+
+# @jit(nopython=True)
+def solve_positions(particles, V):
+    # update the position of each particle from particles according to the velocity acting on it
+    particles_new = particles + dt * V
+    particles_new[:, 1] = np.where((particles_new[:, 1] < 0), 0, particles_new[:, 1]) # if we cross into the floor, reverse the y velocity and apply a damping factor
+    # CANTILEVERED BEAM - set all particles in the first column to have a position of zero
+    '''
+    for particle_num in range(len(particles)):
+        # if the particle is in the first row:
+        if particle_num <= y_num:
+            # keep that particle in place
+            particles_new[particle_num, 0] = particles[particle_num, 0]
+            particles_new[particle_num, 1] = particles[particle_num, 1]
+    '''
+    return particles_new
+
+# Constants:
+x_num = 40 # 10 or 40
+y_num = 40 # 10 or 40
+t_final = 1.0
+dt = 0.005
+sub_increments = 20 # how many timesteps to calculate between each frame
+mesh_spacing = 1
+neutral_distance = mesh_spacing
+height_offset = 8 # 3 for 10x10, 8 for 40x40
+
+# Physical Constants (these alter the particle physics)
+radius_check = neutral_distance*1.35
+spring_constant = 2000 # 1000
+damping_constant = 18 # 10
+damping_constant_global = 0.25 # air resistance
+coef_of_restitution = .1 # 0.5
+gravity = 1
+particle_mass = 0.5
+
+# Create our triangular mesh of points
+particles = np.zeros((x_num*y_num, 2))
+velocities = np.zeros((x_num*y_num, 2))
+forces = np.zeros((x_num*y_num, 2))
+
+# create a mesh of particles, indexes start at point 0,0 and grow up the y-axis, then move to the next column and grow up that one
+for i in range(x_num):
+    for j in range(y_num):
+        particles[i*y_num+j, 0] = i
+        particles[i*y_num+j, 1] = j * np.sqrt(3) / 2 # height of equilateral triangle
+        if j % 2 == 1:  # add an offset to every other point
+            particles[i*y_num+j, 0] += 0.5
+
+# particles[-1, 0] -= mesh_spacing*.5 # offset a particle for chaos
+particles *= mesh_spacing # scale the mesh to the desired spacing
+particles[:,1] += mesh_spacing*height_offset
+centroid = np.mean(particles, axis=0) # find the centroid of the mesh
+
+# rotate the mesh of particles about the origin
+theta = np.radians(30)
+rot_matrix = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])
+particles = np.dot(particles-centroid, rot_matrix) + centroid
+
+# Create a scatterplot of X
+my_dpi = 96
+fig = plt.figure(figsize=(800/my_dpi, 400/my_dpi), dpi=my_dpi)
+ax = fig.add_axes([0, 0, 1, 1], frameon=False)
+# numpy array of colors for each particle in shades from red to blue
+colors = np.zeros((x_num*y_num, 4))
+colors[:, 0] = np.linspace(1, 0, x_num*y_num)
+colors[:, 2] = np.linspace(0, 1, x_num*y_num)
+colors[:, 3] = 1
+plt.axis('equal')   # ensure axes are equally scaled
+plt.plot([-mesh_spacing*10, x_num*mesh_spacing+mesh_spacing*10], [0, 0], color = 'k') # Draw a line representing the ground at y = 0;
+plot = ax.scatter(particles[:, 0], particles[:, 1],
+                  marker = 'o',
+                  s = 20,
+                  color = colors,
+)
+
+time1 = time.time()
+
+
+# @jit(nopython=True)
+def animate(w):
+    global particles, velocities, forces
+    for i in range(sub_increments):
+        # Semi-implicit Euler integration
+#        velocities = velocities + solve_velocities(forces, particles)
+#        particles = solve_positions(particles, velocities)
+#        forces = solve_forces(particles, velocities)
+        # Velocity Verlet integration
+        velocities_half = velocities + solve_velocities(forces, particles, dt/2)
+        particles = solve_positions(particles, velocities_half)
+        forces = solve_forces(particles, velocities)
+        velocities = velocities_half + solve_velocities(forces, particles, dt/2)
+    plot.set_offsets(np.c_[particles[:, 0],particles[:, 1]])
+
+anim = animation.FuncAnimation(fig, animate, frames=int(t_final/dt/2))
+# plt.show()
+solve_time = time.time() - time1
+
+
+time2 = time.time()
+import matplotlib as mpl 
+mpl.rcParams['animation.ffmpeg_path'] = r'C:\ffmpeg_2023_03_30\bin\ffmpeg.exe'
+path = r"c://Users/david/Desktop/NMM/week7/figures/animation.mp4" 
+writervideo = animation.FFMpegWriter(fps=60, extra_args=['-vcodec', 'libx264']) 
+anim.save(path, writer=writervideo, dpi=my_dpi)
+mp4_time = time.time() - time2
+
+print("Solve time: ", solve_time)
+print("MP4 time: ", mp4_time)
\ No newline at end of file