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:

[![enter image description here][1]][1]


  [1]: https://i.sstatic.net/F6Vsq.png