You are not updating the angle in your simulation as the velocity vector is changing. Here is a modified version of your 'animate' function and the result.
def animate(i):
t.append(t[i] + dt)
vx.append(vx[i] + dt * ax[i]) # Update the velocity
vy.append(vy[i] + dt * ay[i])
vel = np.sqrt(vx[-1] ** 2 + vy[-1] ** 2) # magnitude of velocity
cos = vx[-1]/vel
sin = vy[-1]/vel
# Update position
x.append(x[i] + dt * vx[i])
y.append(y[i] + dt * vy[i])
# Calculate updated velocity
drag = area * density * 0.5 * cd * vel ** 2
ax.append(-(drag * cos) / M)
ay.append(-g - (drag * sin / M))
line.set_xdata(x[:-1])
line.set_ydata(y[:-1])
return line,
and the result:
