How do I Decrease Orbit Size in My Python Code?

Question

I'm working on a project where my end goal is for the two orbiting stars to collide; I have some code from my class that I've been utilizing, and while I was successful in getting two objects to orbit one another/orbit the same point, I can't get the radius of the orbit to decrease over time. I'm attempting to add a drag force, but this ejects my stars from their orbits.

G = 6.67e-11 #Actual gravitational constant
    dt = 1e-3
    mS = 1e6
    mp1 = 1e5
    mp2 = 1e5
    x1 = 1 #starting x position
    y1 = 0 #starting y position
    x2 = -1
    y2 = 0
    vx = 0 #starting x-velocity
    vy = 2 #starting y-velocity
    
    for param in ['figure.facecolor', 'axes.facecolor', 'savefig.facecolor']:
        plt.rcParams[param] = 'black' 

    timestep = 500 #timestep for video

    xlist1 = [x1] #x-position of BH 1
    ylist1 = [y1] #y-position of BH 1

    xlist2 = [x2] #x-position of BH 2
    ylist2 = [y2] #y-position of BH 2

    vxlist1 = [vx]
    vylist1 = [vy]

    vxlist2 = [vx]
    vylist2 = [vy]

    n = 10
    
    plt.style.use("cyberpunk")

    for param in ['figure.facecolor', 'axes.facecolor', 'savefig.facecolor']:
        plt.rcParams[param] = 'black' 
        
    tail_length = 100

    for t in range(1, timestep+1):

        for _ in range(n):
            r21 = x1**2 + y1**2
            x1 = xlist1[-1] + vxlist1[-1] * dt #updates x by taking the final value in the xlist array and adding the final value in the velocity (x direction) array multiplied by our time step
            y1 = ylist1[-1] + vylist1[-1] * dt
            vx = vxlist1[-1] - (G*mS*mp1*xlist1[-1]/r21**1.5 + drag force) * dt #gravitational force GM1M2/r^2
            vy = vylist1[-1] - (G*mS*mp1*ylist1[-1]/r21**1.5 + drag force) * dt #update velocity by taking the final value in the vy array and subtracting the gravitational force at the final location logged in the y position array
            xlist1.append(x1)
            ylist1.append(y1)
            vxlist1.append(vx)
            vylist1.append(vy)
        
            r22 = x2**2 + y2**2
            x2 = xlist2[-1] + vxlist2[-1] * dt #updates x by taking the final value in the xlist array and adding the final value in the velocity (x direction) array multiplied by our time step
            y2 = ylist2[-1] + vylist2[-1] * dt
            vx = vxlist2[-1] - (G*mS*mp2*xlist2[-1]/r22**1.5 + drage force) * dt #gravitational force GM1M2/r^2
            vy = vylist2[-1] - (G*mS*mp2*ylist2[-1]/r22**1.5 + drag force) * dt #update velocity by taking the final value in the vy array and subtracting the gravitational force at the final location logged in the y position array
            xlist2.append(x2)
            ylist2.append(y2)
            vxlist2.append(vx)
            vylist2.append(vy)

With this, I've changed velocity, increased and decreased the gravitational constant, and manually tried to decrease the x and y positions by subtracting 0.0001 from each value of the x and y array. At this point, I'm somewhat out of ideas and starting from scratch isn't really an option, so any and all help is very much appreciated :)