diff --git a/kutta.py b/kutta.py
new file mode 100644
index 0000000..d261f39
--- /dev/null
+++ b/kutta.py
@@ -0,0 +1,116 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def RK_n(y0, t, n=4):
+    ODE = lambda t, y: t**3 + y**3
+    #note: t is a vector of time points, y0 is the initial value
+    steps = len(t)
+    y = np.zeros(steps)
+    h = (t[-1] - t[0]) / (steps - 1)
+    if t[0] > 0:
+        raise ValueError("t[0] must be less than or equal to 0")
+    
+    #find the index of the first non-negative time point
+    t0_idx = np.where(t >= 0)[0][0]
+    y[t0_idx] = y0  # set the initial condition at t[0]
+    
+    #forward RK
+    for i in range(t0_idx, steps - 1):
+        #n-th order Runge-Kutta method
+        k = np.zeros(n)
+        k[0] = h * ODE(t[i], y[i])
+        for j in range(1, n):
+            t_j = t[i] + h * (j / n)
+            y_j = y[i] + k[j - 1] / n
+            k[j] = h * ODE(t_j, y_j)
+        y[i + 1] = y[i] + np.sum(k) / n
+        
+    #backward RK
+    for i in range(t0_idx - 1, -1, -1):
+        #n-th order Runge-Kutta method
+        k = np.zeros(n)
+        k[0] = h * ODE(t[i], y[i])
+        for j in range(1, n):
+            t_j = t[i] + h * (j / n)
+            y_j = y[i] + k[j - 1] / n
+            k[j] = h * ODE(t_j, y_j)
+        y[i] = y[i + 1] - np.sum(k) / n
+    return y
+    
+case_num = 3
+if case_num == 1:
+    plt.title("RK4 result with IVPs for large t")
+    initial_values = [-1, 0, 1]
+    t = np.linspace(0, 2, 1000)
+    for i, y0 in enumerate(initial_values):
+        res = RK_n(y0, t, 4)
+        plt.plot(t, res, label=f'y(0)={y0}, n=4')
+    plt.ylim(-5,5)
+    #observation: the solutions diverge as t increases and seem to have an asymptote which depends on y(0)
+    # and converge to a same solution as t decreases.
+    
+elif case_num == 2:
+    plt.title("RK4 result with IVPs for small t")
+    initial_values = [-1, 0, 1]
+    t = np.linspace(-0.5, 0.5, 1000)
+    for i, y0 in enumerate(initial_values):
+        res = RK_n(y0, t, 4)
+        plt.plot(t, res, label=f'y(0)={y0}, n=4')
+    plt.ylim(-5,5)
+
+elif case_num == 3:
+    plt.title("Comparing different order RK methods")
+    order = [1, 2, 4, 5]
+    t = np.linspace(-1, 1.4, 1000)
+    res_5 = RK_n(0, t, 5)
+    for n in order:
+        res = RK_n(0, t, n)
+        if n == 1:
+            label = 'Euler - RK5'
+        elif n == 2:
+            label = 'Improved Euler - RK5'
+        else:
+            label = f'RK{n} - RK5'
+        plt.plot(t, res - res_5, label=label)
+    plt.legend()
+    
+elif case_num == 4:
+    plt.title("RK4 with IVP from 0 to 1")
+    initial_values = np.linspace(0, 1, 5)
+    t = np.linspace(-3, 10, 1000)
+    for i, y0 in enumerate(initial_values):
+        res = RK_n(y0, t, 4)
+        plt.plot(t, res, label=f'y(0)={y0}')
+    plt.ylim(-5,10)
+    plt.legend()
+    
+elif case_num == 5:
+    fig, axs = plt.subplots(2, 2, figsize=(10, 8))
+    plt.title("RK4 with different y-scale in the negatives")
+    scale = [5, 10, 20, 50]
+    t = np.linspace(-10, 1, 1000)
+    res = RK_n(0, t, 4)
+    for i, s in enumerate(scale):
+        ax = axs[i // 2, i % 2]
+        ax.plot(t, res / s, label=f'scale={s}')
+        ax.set_ylim(-1, scale[i])
+        ax.axhline(0, color='black', lw=0.5)
+        ax.axvline(0, color='black', lw=0.5)
+        ax.set_title(f'RK4 for t < {s}')
+    # observation: There is no asymptote in the negative y-scale
+    
+elif case_num == 6:
+    plt.title("RK4 with different step sizes")
+    steps = [5, 10, 100, 1000]
+    for n in steps:
+        t = np.linspace(0, 1, n)
+        res = RK_n(0, t, 4)
+        plt.plot(t, res, label=f'{n} steps')
+    plt.legend()
+    
+
+plt.axhline(0, color='black', lw=0.5)
+plt.axvline(0, color='black', lw=0.5)
+plt.show()
+# plt.legend()
diff --git a/numeric.py b/numeric.py
new file mode 100644
index 0000000..8dfdd03
--- /dev/null
+++ b/numeric.py
@@ -0,0 +1,98 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+threshold = 100
+
+def RK_n(y0, t: np.ndarray, n=4):
+    ODE = lambda t, y: t**3 + y**3
+    #note: t is a vector of time points, y0 is the initial value
+    steps = len(t)
+    y = np.zeros(steps)
+    h = (t[-1] - t[0]) / (steps - 1)
+    if t[0] > 0:
+        raise ValueError("t[0] must be less than or equal to 0")
+    
+    #find the index of the first non-negative time point
+    t0_idx = np.where(t >= 0)[0][0]
+    y[t0_idx] = y0  # set the initial condition at t[0]
+    
+    #forward RK
+    for i in range(t0_idx, steps - 1):
+        #n-th order Runge-Kutta method
+        k = np.zeros(n)
+        k[0] = h * ODE(t[i], y[i])
+        for j in range(1, n):
+            t_j = t[i] + h * (j / n)
+            y_j = y[i] + k[j - 1] / n
+            k[j] = h * ODE(t_j, y_j)
+        y[i + 1] = y[i] + np.sum(k) / n
+    return y
+
+def euler(f, a, b, y0, h):
+    xs = [a]
+    ys = [y0]
+    x = a
+    y = y0
+    while x < b:
+        y += h * f(x, y)
+        x += h
+        if abs(y) > threshold:
+            break
+        xs.append(x)
+        ys.append(y)
+    return xs, ys
+
+def improvedEuler(f, a, b, y0, h):
+    xs = [a]
+    ys = [y0]
+    
+    x = a
+    y = y0
+    while x <= b:
+        slope1 = f(x, y)
+        y_predictor = y + h * slope1
+        slope2 = f(x + h, y_predictor)
+        y = y + (h/2) * (slope1 + slope2)
+        x = x + h
+        if abs(y) > threshold:
+            break
+        xs.append(x)
+        ys.append(y)
+    return xs, ys
+
+lims = [
+    (0, 1.4, 1.8, 1, threshold),
+    (1, 0.4, 0.6, 1, threshold),
+    (-1, 0.4, 0.6, -threshold, -1),
+]
+
+ODE = lambda t, y: t**3 + y**3
+fig, axes = plt.subplots(3, 3)
+
+col_titles = ["Euler's Method", "Improved Euler's", "Runge-Kutta"]
+row_titles = ["y(0) = 0", "y(0) = 1", "y(0) = -1"]
+
+# Add column titles
+for col_idx, col_title in enumerate(col_titles):
+    axes[0][col_idx].set_title(col_title, fontsize=12)
+
+for row_idx, row_title in enumerate(row_titles):
+    axes[row_idx][0].set_ylabel(row_title, fontsize=12, rotation=90, labelpad=20)
+
+for idx, (y0, xm, xM, ym, yM) in enumerate(lims):
+    for h in [0.1, 0.01, 0.001, 0.0001]:
+        Te, Ye = euler(ODE, 0, xM, y0, h)
+        Ti, Yi = improvedEuler(ODE, 0, xM, y0, h)
+        Tr = np.linspace(0, xM, int((xM - 0)/h))
+        Yr = RK_n(y0, Tr)
+        
+        axes[idx][0].plot(Te, Ye, label=f'$\Delta t$={h}', linewidth=1)
+        axes[idx][1].plot(Ti, Yi, label=f'$\Delta t$={h}', linewidth=1)
+        axes[idx][2].plot(Tr, Yr, label=f'$\Delta t$={h}', linewidth=1)
+    for idy in range(3):
+        axes[idx][idy].set_xlim(xm, xM)
+        axes[idx][idy].set_ylim(ym, yM)
+        axes[idx][idy].legend()
+
+fig.suptitle("Numeric Solutions for Each ODE Near the Asymptote")
+plt.show()
diff --git a/picard-poly.py b/picard-poly.py
index 331dd44..1d9acbb 100644
--- a/picard-poly.py
+++ b/picard-poly.py
@@ -188,35 +188,36 @@ def picardIteration(t0, y0, num_iters = 20):
     return phi
 
 lims = [
-    (0, 2, 0, 100),
-    (0, 1, 1, 100),
-    (0, 1, -1, 100)
+    (0, 2, 0, -1, 100),
+    (0, 1, 1, 0, 100),
+    (0, 1, -1, -100, 0)
 ]
 num_iter = 6
 
-t0, tM, y0, yM = lims[1]
-tm, ym = min(t0, tM), min(y0, yM)
+fig, axes = plt.subplots(2, 3)
+for idx, (t0, tM, y0, ym, yM) in enumerate(lims):
+    tm = min(t0, tM)
+    picard_result = picardIteration(t0, y0, num_iter)
+    for idp, poly in enumerate(picard_result):
+        np_poly = np.vectorize(poly)
+        ts = np.linspace(t0, tM, 1000)
+        ys = np_poly(ts)
+        axes[0][idx].plot(ts, ys, label=f'#{idp} iteration')
+    axes[0][idx].legend()
+    axes[0][idx].set_xlim(tm, tM)
+    axes[0][idx].set_ylim(ym, yM)
+    
+    #last iter coeffs
+    axes[1][idx].scatter(list(range(poly.degree + 1)), poly.coefficients, s=0.5)
+    axes[1][idx].set_yscale('log')
 
-picard_result = picardIteration(t0, y0, num_iter)
-'''
-for idx, poly in enumerate(picard_result):
-    np_poly = np.vectorize(poly)
-    ts = np.linspace(t0, tM, 1000)
-    ys = np_poly(ts)
-    plt.plot(ts, ys)
-    print(idx, poly.degree, poly)
-'''
-poly = picard_result[-1]
-#np_poly = np.vectorize(poly)
-#ts = np.linspace(t0, tM, 1000)
-#ys = np_poly(ts)
-#plt.plot(ts, ys)
-#print(poly.degree, poly)
-print(poly.degree)
-plt.scatter(list(range(poly.degree + 1)), poly.coefficients, s=0.5)
-plt.yscale('log')
-#plt.xlim(t0, tM)
-#plt.ylim(ym, yM)
-plt.title(f"coefficients of $\\sum c_n t^n$ at y(0) = {y0} of iterations {num_iter}")
-plt.ylabel("$c_n$")
-plt.show()
\ No newline at end of file
+col_titles = ["y(0) = 0", "y(0) = 1", "y(0) = -1"]
+row_titles = ["plots of solutions", "coefficients for $t^n$"]
+
+# Add column titles
+for col_idx, col_title in enumerate(col_titles):
+    axes[0][col_idx].set_title(col_title, fontsize=12)
+
+for row_idx, row_title in enumerate(row_titles):
+    axes[row_idx][0].set_ylabel(row_title, fontsize=12, rotation=90, labelpad=20)
+plt.show()