I'll give you an intuitive way to think about it, then a script to play around with if you like.
When the spacecraft is at periapsis, it is there, at periapsis. A bit of drag isn't going to change its position by much, but it will change its velocity.
The new velocity determines how high it can reach at apoapsis. If it has slowed down, it will climb to a lower apoapsis.
However, it's in an elliptical orbit, and there's no drag up there, so it's going to return to nearly the same periapsis whence it came.
caveat: That assumes there is no lift. The term "lift" stands for any aerodynamic force not in the direction of the relative velocity of the spacecraft with respect to the atmosphere. If the spacecraft is not spherical or sufficiently tumbling, or if it has wings and is space-shuttle-shaped and Clint Eastwood is in the "pilot's seat", then all bets are off.
Here's a python script for a simple approximate simulation that demonstrates this. It is not intended to be precise, but it's physical enough to show that this is what tends to happen.
Eccentric orbits "touching the atmosphere" at periapsis will tend to first circularize before they burn up.
And this has at least one important implication; it won't necessarily re-enter and burn up near the original periapsis. It can happen anywhere. Round and round it goes; where it stops, nobody knows.
x, v = X.reshape(2, -1)
alt = np.sqrt((x**2).sum()) - re
rho = rho0 * np.exp(-alt/hscale)
Fdrag = -0.5 * rho * CD * area * v * np.sqrt((v**2).sum())
def deriv(X, t):
x, v = X.reshape(2, -1)
acc_g = -GMe * x * ((x**2).sum())**-1.5
acc_d = acc_drag(X)
return np.hstack((v, acc_d + acc_g))
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint
halfpi, pi, twopi = [f*np.pi for f in [0.5, 1.0, 2.0]]
GMe = 3.98600418E+14 # m^3 s^-2
mass = 5000. # kg
area = 5. # m^2
CD = 1.0
hscale = 7200. # meters (ROUGHLY scale height fudged for 100km)
re = 6378000. # meters (equatorial radius)
rho0 = 1.225 # kg/m^3
r_peri = re + 90000. # meters
v_peri = 10000. # m/s
X0 = np.hstack(((r_peri, 0, 0), (0, v_peri, 0)))
time = np.arange(0, 1000000, 100) # sec
answer, info = ODEint(deriv, X0, time, full_output=True)
theta = np.linspace(0, twopi, 360)
xe, ye = [re*f(theta) for f in [np.cos, np.sin]]
alti = np.sqrt((answer[:,:3]**2).sum(axis=1)) - re
plt.subplot(2, 1, 1)
x, y = answer.T[:2]
plt.plot(, , 'ok')
plt.plot(xe, ye, '-k', linewidth=1.5)
plt.subplot(2, 1, 2)