Thrust = m(0)*a
OK for when drag is low, like at the very beginning or above something like 30 km, as can be seen from your dynamic pressure graph.
So for some time t you'd need to add drag:
Thrust = m(t)*a + (1/2)*rho(t)*v(t)^2*C_d*A. The obvious problem here is that rho depends on altitude (or equivalently on time), and there's speed again, so much more complicated. I'm not sure if you're doing this? Otherwise you end up with artificially lower thrust when drag is significant.

m_dot = (Thrust - Delta_pA)/(ISPg0)
Here you also need to know how pressure depends on altitude (and hence time) and also what's the nozzle exit area. There's also some dependence of Isp on altitude.

But anyway, ignoring nuances, the bottom line is, for accurate thrust data you definitely need to factor in drag around maxq.