diff --git a/README.md b/README.md
index 8e9c04d..4ddc63c 100644
--- a/README.md
+++ b/README.md
@@ -93,3 +93,47 @@ function y = prob3ptC(prob3_ode, dt, y0,tlimits)
 end
 ```
 ![](assets/README-fcbec9a1.png)
+
+# Problem 4
+## Part A
+The analytical solution to v'=g-cd/m v*v is
+
+## Part B
+```
+cd = 0.25;
+m = 60;
+g = 9.81;
+dt = 0.1;
+t = [0:dt:12]';
+x = zeros(length(t),2);
+x(1,:) = [1000, 0];
+for i = 2:length(t)
+    dx = prob4_ode(t,x(i-1,:));
+    x(i,:) = x(i-1,:)+dt*dx;
+end
+```
+
+## Part C
+```
+function x = prob4ptC(prob4_ode, dt, x0,tlimits)
+    t = [tlimits(1):dt:tlimits(2)];
+    x = zeros(length(t),length(x0));
+    x(1,:) = x0;
+    for i = 2:length(t)
+        dx = prob4_ode(t(i-1),x(i-1,:));
+        x(i,:) = x(i-1,:)+dt*dx;
+        dxc = (dx + prob4_ode(t(i),x(i,:)))/2;
+        x(i,:) = x(i-1,:)+dt*dxc;
+        n = 1;
+        while (1)
+            n = n+1;
+            yold = x(i,:);
+            dxc = (dx+prob4_ode(t(i),x(i,:)))/2;
+            x(i,:) = x(i-1,:)+dt*dxc;
+            ea = abs(x(i,:)- yold)./x(i,:)*100;
+            if max(ea) < 0.0001 || n>100; break, end
+        end
+    end
+end
+```
+![](assets/README-ac6c7a25.png)
diff --git a/assets/README-ac6c7a25.png b/assets/README-ac6c7a25.png
new file mode 100644
index 0000000..d5d0d51
Binary files /dev/null and b/assets/README-ac6c7a25.png differ
diff --git a/prob4_ode.m b/prob4_ode.m
new file mode 100644
index 0000000..926a432
--- /dev/null
+++ b/prob4_ode.m
@@ -0,0 +1,8 @@
+function dx = ode4(t,x)
+g = 9.81;
+cd = 0.25;
+m = 60;
+dx = zeros(size(x));
+dx(1) = x(2); 
+dx(2) = g -(cd/m)*x(2)^2;
+end 
\ No newline at end of file
diff --git a/prob4plot.m b/prob4plot.m
new file mode 100644
index 0000000..f41d964
--- /dev/null
+++ b/prob4plot.m
@@ -0,0 +1,20 @@
+set(0, 'defaultAxesFontSize', 16)
+set(0,'defaultTextFontSize',14)
+set(0,'defaultLineLineWidth',3)
+cd = 0.25;
+m = 60;
+g = 9.81;
+dt = 0.1;
+t = [0:dt:12]';
+x = zeros(length(t),2);
+x(1,:) = [1000, 0];
+for i = 2:length(t)
+    dx = prob4_ode(t,x(i-1,:));
+    x(i,:) = x(i-1,:)+dt*dx;
+end
+c = figure(3);
+x_heun = prob4ptC (@(t,x) prob4_ode(t,x), 0.1, [1000, 0], [0 12]);
+plot(t, m/cd*log(cosh(sqrt((g*cd)/m)*t))+1000,t,x(:,1),t,x_heun(:,1)); 
+legend('Analytical','Euler','Heuns','Location','Northeast')
+xlabel('Time')
+ylabel('Position')
\ No newline at end of file
diff --git a/prob4ptC.m b/prob4ptC.m
new file mode 100644
index 0000000..f1dff21
--- /dev/null
+++ b/prob4ptC.m
@@ -0,0 +1,20 @@
+function x = prob4ptC(prob4_ode, dt, x0,tlimits)
+    t = [tlimits(1):dt:tlimits(2)];
+    x = zeros(length(t),length(x0));
+    x(1,:) = x0;
+    for i = 2:length(t)
+        dx = prob4_ode(t(i-1),x(i-1,:));
+        x(i,:) = x(i-1,:)+dt*dx;
+        dxc = (dx + prob4_ode(t(i),x(i,:)))/2;
+        x(i,:) = x(i-1,:)+dt*dxc;
+        n = 1;
+        while (1) 
+            n = n+1;
+            yold = x(i,:);
+            dxc = (dx+prob4_ode(t(i),x(i,:)))/2;
+            x(i,:) = x(i-1,:)+dt*dxc;
+            ea = abs(x(i,:)- yold)./x(i,:)*100; 
+            if max(ea) < 0.0001 || n>100; break, end
+        end
+    end
+end 
\ No newline at end of file