import numpy as np
import matplotlib.pyplot as plt

# Problem1
h = [0.001, 0.005, 0.01, 0.05, 0.1]
num = 5 / h[0]
y = 1
t = 0
dy_dt = -1
tPlot = np.arange(0, 5, h[0])
yPlot = []
for i in range(int(num)):
    yPlot.append(y)
    t += h[0]
    y += dy_dt * h[0]
    dy_dt = -2*y+2-np.exp(-4*t)
    if i % 1000 == 999:
        print(t, y)
plt.plot(tPlot, yPlot, label="Euler's method")
yExact = (np.exp(-4*tPlot))/2-(np.exp(-2*tPlot))/2+1
plt.plot(tPlot, yExact, label="exact solution")
plt.xlabel("t")
plt.ylabel("y")
plt.ylim(0.8, 1.05)
plt.legend()
plt.show()

plt.plot(tPlot, yPlot, label="0.001")
for i in range(4):
    num = 5 / h[i+1]
    y = 1
    t = 0
    dy_dt = -1
    tPlot = np.arange(0, 5, h[i+1])
    yPlot = []
    for j in range(int(num)):
        yPlot.append(y)
        t += h[i+1]
        y += dy_dt * h[i+1]
        dy_dt = -2*y+2-np.exp(-4*t)
    plt.plot(tPlot, yPlot, label=h[i+1])
plt.xlabel("t")
plt.ylabel("y")
plt.ylim(0.8, 1.05)
plt.legend()
plt.show()