diff --git a/week4/README.MD b/week4/README.MD
new file mode 100644
index 0000000000000000000000000000000000000000..be5aeba2ed236a61a7e5a712e954b1f227d9fdf4
--- /dev/null
+++ b/week4/README.MD
@@ -0,0 +1,9 @@
+# Week 4 - Ordinary Differential Equations
+
+From 6.2, error at different step sizes for Euler integration (both are diverging) and Runge-Kutta with much smaller error.
+
+
+
+From 6.3, adaptive Runge-Kutta stepper's step size increasing to meet a desired error target. Could have used a more aggressive minimum error term to see the steps oscillate with the derivative instead of just meeting the requirements of the peak derivative.
+
+
diff --git a/week4/week4_2.py b/week4/week4_2.py
new file mode 100644
index 0000000000000000000000000000000000000000..9db6d6257a6070659631eea630414efd20bc4626
--- /dev/null
+++ b/week4/week4_2.py
@@ -0,0 +1,64 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+'''
+6.2)
+For a simple harmonic oscillator y'' + y = 0, with initial conditions y(0) = 1, y'(0) = 0,
+find y(t) from t = 0 to 100π. Use an Euler method and a fixed-step fourth-order Runge-Kutta method.
+For each method check how the average local error, and the error in the final value and slope, depend on the step size
+'''
+
+# Euler method
+def euler_integrate(step):
+ t = np.arange(0, 100*np.pi, step)
+ increments = np.arange(1, len(t), 1) # we will need something indexed from 1 to final time_step to not interfere with the initial conditions
+ y = np.zeros(len(t))
+ y[0] = 1
+ dy = np.zeros(len(t))
+ dy[0] = 0
+ for time in increments:
+ y[time] = y[time-1] + step*dy[time-1]
+ dy[time] = dy[time-1] - step*(y[time-1])
+ # dy[time] = dy[time-1] - step*(y[time]) # reverse euler? Just use the expected new y
+ return y, dy, t
+
+# Runge-Kutta method
+def rk4_integrate(step):
+ t = np.arange(0, 100*np.pi, step)
+ increments = np.arange(1, len(t), 1) # we will need something indexed from 1 to final time_step to not interfere with the initial conditions
+ y = np.zeros(len(t))
+ y[0] = 1
+ dy = np.zeros(len(t))
+ dy[0] = 0
+ for time in increments:
+ k1, k2, k3, k4 = [0, 0], [0, 0], [0, 0], [0, 0] # make position and velocity vectors
+
+ k1[0] = step * dy[time-1] # normal euler step
+ k1[1] = - step * y[time-1]
+ k2[0] = step * (dy[time-1] + k1[1]*0.5) # half euler step in the direction of k1
+ k2[1] = -step * (y[time-1] + k1[0]*0.5)
+ k3[0] = step * (dy[time-1] + k2[1]*0.5) # half euler step in the direction of k2
+ k3[1] = -step * (y[time-1] + k2[0]*0.5)
+ k4[0] = step * (dy[time-1] + k3[1]) # full euler step in the direction of k3
+ k4[1] = -step * (y[time-1] + k3[0])
+
+ y[time] = y[time-1] + (k1[0] + 2* k2[0] + 2 * k3[0] + k4[0])/6 # divide by 6 (not 4) because we assigned 2x value to k2 and k3
+ dy[time] = dy[time-1] + (k1[1] + 2* k2[1] + 2 * k3[1] + k4[1])/6
+ return y, dy, t
+
+step1 = 0.001
+step2 = 0.01
+
+# y1, dy1, t1 = euler_integrate(step1)
+# y2, dy2, t2 = euler_integrate(step2)
+
+y1, dy1, t1 = rk4_integrate(step1)
+y2, dy2, t2 = rk4_integrate(step2)
+
+fig, ax = plt.subplots(2)
+ax[0].plot(t1, y1)
+ax[0].title.set_text('Times step = ' + str(step1) + '\n Final y: ' + str(round(y1[-1],3)) + '\n Final dy: ' + str(round(dy1[-1],3)))
+ax[1].plot(t2, y2)
+ax[1].title.set_text('Times step = ' + str(step2) + '\n Final y: ' + str(round(y2[-1],3)) + '\n Final dy: ' + str(round(dy2[-1],3)))
+plt.tight_layout()
+plt.show()
\ No newline at end of file
diff --git a/week4/week4_3.py b/week4/week4_3.py
new file mode 100644
index 0000000000000000000000000000000000000000..c79fa9fd0ce2263144911bc64acd917f0d933d55
--- /dev/null
+++ b/week4/week4_3.py
@@ -0,0 +1,73 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+'''
+6.3)
+Write a fourth-order Runge Kutta adaptive stepper for the preceding problem,
+and check how the average step size that it finds depends on the desired local error.
+'''
+
+# Adaptive Runge-Kutta method (rewritten to be a step call)
+def rk4_step(step, y, dy):
+ k1, k2, k3, k4 = [0, 0], [0, 0], [0, 0], [0, 0] # make position and velocity vector
+ k1[0] = step * dy # normal euler step for velocity
+ k1[1] = - step * y # normal euler step for position
+ k2[0] = step * (dy + k1[1]*0.5) # half euler step in the direction of k1
+ k2[1] = -step * (y + k1[0]*0.5)
+ k3[0] = step * (dy + k2[1]*0.5) # half euler step in the direction of k2
+ k3[1] = -step * (y + k2[0]*0.5)
+ k4[0] = step * (dy + k3[1]) # full euler step in the direction of k3
+ k4[1] = -step * (y + k3[0])
+ y_new = y + (k1[0] + 2*k2[0] + 2*k3[0] + k4[0])/6 # divide by 6 (not 4) because we assigned 2x value to k2 and k3
+ dy_new = dy + (k1[1] + 2*k2[1] + 2*k3[1] + k4[1])/6
+ return y_new, dy_new
+
+def rk4_adaptive(step_0, step_nudge, error_min, error_max):
+ ''' In this case we don't know how many timesteps we will take, so it's best to
+ start with lists and convert to numpy arrays at the end...
+ '''
+ t_final = 100*np.pi
+ t = [0]
+ y = [1]
+ dy = [0]
+ steps = [0.1]
+
+ step = steps[0]
+ time = step
+ count = 1
+
+ while time < t_final:
+ y_full, dy_full = rk4_step(step, y[count-1], dy[count-1])
+ y_half1, dy_half1 = rk4_step(step*0.5, y[count-1], dy[count-1])
+ y_half2, dy_half2 = rk4_step(step*0.5, y_half1, dy_half1)
+
+ error = np.abs(y_full - y_half2) # error is how closely the full and half steps agree
+ # if we are too high or too low, nudge the step size and try again
+ if error > error_max:
+ step *= step_nudge
+ elif error < error_min:
+ step *= 1 - (step_nudge-1)
+ # otherwise go ahead and append and increment time and count
+ else:
+ y.append(y_full)
+ dy.append(dy_full)
+ steps.append(step)
+ time += step
+ t.append(time)
+ count += 1
+ return t, y, dy, steps
+
+step_0 = 0.1
+nudge = 0.99 # increase or decrease the step by this percentage to nudge the step size
+error_min = 1e-6 # easier to choose if theses are informed by something physical
+error_max = 1e-4
+
+t, y, dy, steps = rk4_adaptive(step_0, nudge, error_min, error_max)
+
+fig, ax = plt.subplots(2)
+ax[0].plot(t, y)
+ax[0].title.set_text('\n Max Error: ' + str(round(error_max, 10)) + '\n Min Error: ' + str(round(error_min, 10)))
+ax[1].plot(t, steps)
+ax[1].title.set_text('Total number of steps: ' + str(len(t)) + '\n Average step size: ' + str(round(np.mean(steps), 6)))
+plt.tight_layout()
+plt.show()
diff --git a/week4/week4_4.py b/week4/week4_4.py
new file mode 100644
index 0000000000000000000000000000000000000000..41b64b38ffcebfa14179fc4ebbb1817e8c0f74d0
--- /dev/null
+++ b/week4/week4_4.py
@@ -0,0 +1,65 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+'''
+6.4)
+Numerically solve the differential equation found in Problem 4.3:
+l*ddθ + (g + ddz)*sinθ = 0
+Take the motion of the platform to be periodic, and interactively explore the dynamics of
+the pendulum as a function of the amplitude and frequency of the excitation.
+
+Consider a pendulum with a mass m connected by a massless rod of length l to a moveable platform at height z (Figure 6.4)
+'''
+step = 0.001
+t_final = 10*np.pi
+t = np.arange(0, t_final, step)
+increments = np.arange(1, len(t), 1)
+time = 0
+
+# Platform Parameters (periodic motion and very differentiable)
+A = 1 # amplitude of the platform's motion [m]
+platform_omega = 1 # angular frequency of the platform's motion [rad/s]
+zs = A*np.sin(platform_omega*t) # all platform positions as a function of time
+ddzs = A*platform_omega**2*np.cos(platform_omega*t) # all platform accelerations as a function of time
+
+# Pendulum Parameters
+thetas = np.zeros(len(t))
+d_thetas = np.zeros(len(t))
+dd_thetas = np.zeros(len(t))
+l = 1 # length of pendulum
+g = 9.81 # acceleration due to gravity
+theta_0 = 0 # initial angle of the pendulum [rad]
+d_theta_0 = 1 # initial angular velocity of the pendulum [rad/s]
+
+def d_theta(theta, omega, t):
+ return omega
+
+def dd_theta(theta, omega, t):
+ solution = -np.sin(theta)/l * (g - ddzs[int(t/step)])
+ print(solution)
+ return solution
+
+def rk4_step(prev_values, dif_eqs, step, time):
+ variables = np.array(prev_values)
+ k1, k2, k3, k4 = [0, 0], [0, 0], [0, 0], [0, 0] # make position and velocity vector
+ ddy = dif_eqs[1](prev_values[0], prev_values[1], time)
+ dy = dif_eqs[0](prev_values[0], prev_values[1], time)
+ print(ddy)
+ k1[0] = step * ddy # normal euler step for velocity
+ k1[1] = - step * dy # normal euler step for position
+ k2[0] = step * (ddy + k1[1]*0.5) # half euler step in the direction of k1
+ k2[1] = -step * (dy + k1[0]*0.5)
+ k3[0] = step * (ddy + k2[1]*0.5) # half euler step in the direction of k2
+ k3[1] = -step * (dy + k2[0]*0.5)
+ k4[0] = step * (ddy + k3[1]) # full euler step in the direction of k3
+ k4[1] = -step * (dy + k3[0])
+ variables = variables + np.array(k1)/6 + np.array(k2)/3 + np.array(k3)/3 + np.array(k4)/6
+ return variables
+
+for increment in increments:
+ time = t[increment]
+ thetas[increment], d_theta[increment] = rk4_step([thetas[increment-1], d_thetas[increment-1]], [d_theta, dd_theta], step, time)
+
+# let's plot the chaotic theta over the periodic z
+plt.plot(t, zs)
+plt.plot(t, thetas)
\ No newline at end of file