import matplotlib.pyplot as plt
import numpy as np

dt = 0.003
time = np.arange(0, 100*np.pi, dt)
n_steps = len(time)

#analytical solution
y_correct = np.cos(time)

#initialize y
y = np.zeros(n_steps)
#initialize 1st derivative of y
y_ = np.zeros(n_steps)
#initialize 2nd derivative of y
y__ = np.zeros(n_steps)

error = np.zeros(n_steps)

#initial conditions
y[0] = 1
y_[0] = 0
y[1] = 1
y_[1] = 0
error[0] = 0
error[1] = 0

for i in range(n_steps-2):
    y__[i] = -y[i]
    y_[i+1] = y_[i] + dt*y__[i]
    y[i+2] = y[i+1] + dt*y_[i+1]
    error[i+2] = abs(y[i+2] - y_correct[i+2])

plt.plot(time, error , 'r')
plt.plot(time, y_correct, 'black', linewidth=1)
plt.show()