diff --git a/README.md b/README.md
index c997aa6..8e9c04d 100644
--- a/README.md
+++ b/README.md
@@ -16,9 +16,10 @@ for i=2:length(x)
     dy=prob2_ode(x,y_n(i-1));
     y_n(i)=y_n(i-1)+ dt(2)*dy;
 end  
-plot(x,exp(-x),'k-',x,y_n,'o');
+plot(x,exp(-x),x,y_n)
+legend('Analytical','Euler','Location','Northeast')
 ```
-![](assets/README-66741291.png)
+![](assets/README-be88b588.png)
 
 ## Part C
 
@@ -46,7 +47,7 @@ y(1)=y0;
     end
 end
 ```
-![](assets/README-23f54e97.png)
+![](assets/README-8b7eb1c9.png)
 
 # Problem 3
 ## Part A
@@ -64,6 +65,31 @@ for i=2:length(x)
         dy=prob3_ode(x,y_n(i-1,:));
         y_n(i,:)=y_n(i-1,:)+ dt(2)*dy;
 end
-plot(x,cos(3*x),'k-',x,y_n(:,1),'o');
+plot(x,cos(3*x),x,y_n(:,1));
 ```
-![](assets/README-4e4bb548.png)
+![](assets/README-bec1ab77.png)
+
+## Part C
+```
+function y = prob3ptC(prob3_ode, dt, y0,tlimits)
+    t = [tlimits(1):dt:tlimits(2)];
+    y = zeros(length(t),length(y0));
+    y(1,:) = y0;
+    for i = 2:length(t)
+        dy = prob3_ode(t(i-1),y(i-1,:));
+        y(i,:) = y(i-1,:)+dt*dy;
+        dyc = (dy + prob3_ode(t(i),y(i,:)))/2;
+        y(i,:) = y(i-1,:)+dt*dyc;
+        n = 1;
+        while (1)
+            n = n+1;
+            yold = y(i,:);
+            dyc = (dy+prob3_ode(t(i),y(i,:)))/2;
+            y(i,:) = y(i-1,:)+dt*dyc;
+            ea = abs(y(i,:)- yold)./y(i,:)*100;
+            if max(ea) < 0.0001 || n>100; break, end
+        end
+    end
+end
+```
+![](assets/README-fcbec9a1.png)
diff --git a/assets/README-23f54e97.png b/assets/README-23f54e97.png
deleted file mode 100644
index a16cda4..0000000
Binary files a/assets/README-23f54e97.png and /dev/null differ
diff --git a/assets/README-4e4bb548.png b/assets/README-4e4bb548.png
deleted file mode 100644
index 5dccf9e..0000000
Binary files a/assets/README-4e4bb548.png and /dev/null differ
diff --git a/assets/README-66741291.png b/assets/README-66741291.png
deleted file mode 100644
index 5993ace..0000000
Binary files a/assets/README-66741291.png and /dev/null differ
diff --git a/assets/README-8b7eb1c9.png b/assets/README-8b7eb1c9.png
new file mode 100644
index 0000000..39646e5
Binary files /dev/null and b/assets/README-8b7eb1c9.png differ
diff --git a/assets/README-be88b588.png b/assets/README-be88b588.png
new file mode 100644
index 0000000..a6ebdbd
Binary files /dev/null and b/assets/README-be88b588.png differ
diff --git a/assets/README-bec1ab77.png b/assets/README-bec1ab77.png
new file mode 100644
index 0000000..98d8b73
Binary files /dev/null and b/assets/README-bec1ab77.png differ
diff --git a/assets/README-fcbec9a1.png b/assets/README-fcbec9a1.png
new file mode 100644
index 0000000..5651afb
Binary files /dev/null and b/assets/README-fcbec9a1.png differ
diff --git a/prob2plot.m b/prob2plot.m
index be9521c..4cd4ea4 100644
--- a/prob2plot.m
+++ b/prob2plot.m
@@ -11,7 +11,7 @@ for i = 2:length(t)
 end
 a = figure(1);
 y_heun = Prob2ptC (@(t,y) prob2_ode(t,y), 0.1, 1, [0 3]);
-plot(t,exp(-t),'k-',t,y,'Xg',t,y_heun,'Or'); 
+plot(t,exp(-t),t,y,t,y_heun); 
 legend('Analytical','Euler','Heun','Location','Northeast')
 xlabel('Time')
 ylabel('Position')
\ No newline at end of file
diff --git a/prob2ptB.asv b/prob2ptB.asv
new file mode 100644
index 0000000..43a380e
--- /dev/null
+++ b/prob2ptB.asv
@@ -0,0 +1,12 @@
+dt=[0.5 .1 .001];
+figure();
+hold on
+     x=[0:dt(2):3]';
+     y_n=zeros(size(x));
+    y_n(1)=1;
+for i=2:length(x)
+    dy=prob2_ode(x,y_n(i-1));
+    y_n(i)=y_n(i-1)+ dt(2)*dy;
+end  
+plot(x,exp(-x),x,y_n)
+legend('Analytical','Euler','Location','Northeast')
\ No newline at end of file
diff --git a/prob2ptB.m b/prob2ptB.m
index 4f336d6..43a380e 100644
--- a/prob2ptB.m
+++ b/prob2ptB.m
@@ -8,4 +8,5 @@ for i=2:length(x)
     dy=prob2_ode(x,y_n(i-1));
     y_n(i)=y_n(i-1)+ dt(2)*dy;
 end  
-plot(x,exp(-x),'k-',x,y_n,'o');
\ No newline at end of file
+plot(x,exp(-x),x,y_n)
+legend('Analytical','Euler','Location','Northeast')
\ No newline at end of file
diff --git a/prob3plot.asv b/prob3plot.asv
new file mode 100644
index 0000000..ca98161
--- /dev/null
+++ b/prob3plot.asv
@@ -0,0 +1,17 @@
+set(0, 'defaultAxesFontSize', 16)
+set(0,'defaultTextFontSize',14)
+set(0,'defaultLineLineWidth',1)
+dt = 0.01;
+t = [0:dt:3];
+y = zeros(length(t),2);
+y(1,:) = [1,0];
+for i = 2:length(t)
+    dy = prob3_ode(t,y(i-1,:));
+    y(i,:) = y(i-1,:)+dt*dy;
+end
+b = figure(2);
+y_heun = prob3ptC (@(t,y) prob3_ode(t,y),0.01,[1,0],[0 3]);
+plot(t,cos(3*t),t, y(:,1),t,y_heun(:,1))
+legend('Analytical','Euler','Heun','Location','Northeast')
+xlabel('Time')
+ylabel('Position')
diff --git a/prob3plot.m b/prob3plot.m
new file mode 100644
index 0000000..ca98161
--- /dev/null
+++ b/prob3plot.m
@@ -0,0 +1,17 @@
+set(0, 'defaultAxesFontSize', 16)
+set(0,'defaultTextFontSize',14)
+set(0,'defaultLineLineWidth',1)
+dt = 0.01;
+t = [0:dt:3];
+y = zeros(length(t),2);
+y(1,:) = [1,0];
+for i = 2:length(t)
+    dy = prob3_ode(t,y(i-1,:));
+    y(i,:) = y(i-1,:)+dt*dy;
+end
+b = figure(2);
+y_heun = prob3ptC (@(t,y) prob3_ode(t,y),0.01,[1,0],[0 3]);
+plot(t,cos(3*t),t, y(:,1),t,y_heun(:,1))
+legend('Analytical','Euler','Heun','Location','Northeast')
+xlabel('Time')
+ylabel('Position')
diff --git a/prob3ptB.m b/prob3ptB.m
index ee1e635..bdedc89 100644
--- a/prob3ptB.m
+++ b/prob3ptB.m
@@ -8,4 +8,5 @@ for i=2:length(x)
         dy=prob3_ode(x,y_n(i-1,:));
         y_n(i,:)=y_n(i-1,:)+ dt(2)*dy;
 end
-plot(x,cos(3*x),'k-',x,y_n(:,1),'o');
\ No newline at end of file
+plot(x,cos(3*x),x,y_n(:,1))
+legend('Analytical','Euler','Location','Northeast')
\ No newline at end of file
diff --git a/prob3ptC.m b/prob3ptC.m
new file mode 100644
index 0000000..201174c
--- /dev/null
+++ b/prob3ptC.m
@@ -0,0 +1,20 @@
+function y = prob3ptC(prob3_ode, dt, y0,tlimits)
+    t = [tlimits(1):dt:tlimits(2)];
+    y = zeros(length(t),length(y0));
+    y(1,:) = y0;
+    for i = 2:length(t)
+        dy = prob3_ode(t(i-1),y(i-1,:));
+        y(i,:) = y(i-1,:)+dt*dy;
+        dyc = (dy + prob3_ode(t(i),y(i,:)))/2;
+        y(i,:) = y(i-1,:)+dt*dyc;
+        n = 1;
+        while (1) 
+            n = n+1;
+            yold = y(i,:);
+            dyc = (dy+prob3_ode(t(i),y(i,:)))/2;
+            y(i,:) = y(i-1,:)+dt*dyc;
+            ea = abs(y(i,:)- yold)./y(i,:)*100; 
+            if max(ea) < 0.0001 || n>100; break, end
+        end
+    end
+end
\ No newline at end of file