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:

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 # note: should deal with vel=0 properly
    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:

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:

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: