Wikipedia's Drag equation is
$$F_D = \frac{1}{2} \rho v^2 C_D A$$
shows drag's $\frac{1}{2} \rho v^2$ dependence you mention, as does @MarkAdler's answer about dynamic pressure or "Q":
Max Q is simply the maximum of the dynamic pressure of the external flow, ${1\over 2}\rho v^2$. It has nothing to do with the vehicle, except for the vehicle's speed relative to the undisturbed fluid.
The answer to your question is in their answer. Q's definition does not relate to the body or its aerodynamic properties at all, whereas drag is all about aerodynamic properties!
It's the drag coefficient or $C_D$ that varies, the coefficient of drag.
What is $C_D$ exactly? It's really just a place to put all of the variability in the drag of a real-world object, namely, it's the ratio of the real world (measured or numerically simulated) to $\frac{1}{2} \rho v^2$:
$$C_D = \frac{F_D \ \text{ measured or simulated}}{\frac{1}{2} \rho v^2 A}.$$
If $C_D$ were constant then both $F_D$ and $Q$ would be maximum when $\frac{1}{2} \rho v^2$ reached maximum, but due to the realities of aerodynamics $C_D$ varies quite a lot.
Here's an example of Saturn-V drag coefficient archived from braeunig.us/apollo/saturnV from @RussellBorogove's answer found in this answer:

More generally:

Drag Coefficients as a function of Mach number from source and Robert A. Braeunig's excellent main site
Per request, here's how to simulate this with Python.
Drag and Q look identical because c_D
is a constant. All you have to do is put some functional form for it in the deriv()
function and you can have them reach maxima at different times.
Right now I have the thrust suddenly cut to zero at t_burn
. I've set rtol = 1E-11
to minimize the impact on that non-smooth function. It's better if you ramp thrust down smoothly over say 1 second, and I've shown how to do that in this answer.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
def deriv(t, X):
h, v = X
m = m_f + (m_i - m_f) * np.maximum((t_burn - t) / t_burn, 0)
acc_thrust = (T_avg / m) * (t <= t_burn)
# ugly non-smooth, better use use https://space.stackexchange.com/a/44542/12102
acc_grav = -g
rho = rho_0 * np.exp(-h / H)
F_drag = 0.5 * rho * v**2 * c_D * A
acc_drag = -F_drag / m
acc = acc_thrust + acc_drag + acc_grav
return np.hstack((v, acc))
# from https://space.stackexchange.com/q/51355/12102
g = 9.81 # gravitational acceleration [m/s^2]
rho_0 = 1.225 # air density at sea level [kg/m^3]
H = 8800 # scale height [m]
c_D = 0.54 # drag coefficient for the rocket [-]
A = 0.0103 # rocket body frontal area [m^2]
m_i = 19.1 # “wet” mass of the rocket [kg]
m_f = 10.604 # “dry” mass of the rocket [kg]
T_avg = 2500 # average rocket engine thrust [N]
t_burn = 6.09 # burn time [s]
# Saving some parameters to optimize run time
m = (m_i + m_f)/2 # average mass of the rocket
k = 0.5*c_D*A*rho_0 # constant factor in the drag function
X0 = np.zeros(2)
times = np.linspace(0, 10, 1001)
t_span = times[[0, -1]] # first and last time
answer_drag = solve_ivp(deriv, t_span, X0, t_eval=times, rtol=1E-11)
h, v = answer_drag.y
# Recalculate these at times
mass = m_f + (m_i - m_f) * np.maximum((t_burn - times) / t_burn, 0)
rho = rho_0 * np.exp(-h / H)
F_drag = 0.5 * rho * v**2 * c_D * A
acc_drag = F_drag / mass
Q = 0.5 * rho * v**2
fig, (ax1, ax2, ax3, ax4, ax5, ax6) = plt.subplots(6)
ax1.plot(times, 0.001 * h)
ax1.set_ylabel('height (km)')
ax2.plot(times, v)
ax2.set_ylabel('speed (m/s)')
ax3.plot(times, mass)
ax3.set_ylim(0, None)
ax3.set_ylabel('mass (kg)')
ax4.plot(times, rho)
ax4.set_ylabel('density (kg/m^3)')
ax5.plot(times, F_drag)
ax5.set_ylabel('F_drag (N)')
ax6.plot(times, Q)
ax6.set_ylabel('Q (Pa (kg / m s^2))')
ax6.set_xlabel('time (sec)')
plt.show()