how should I notate my bounds of integration?
I guess you mean for this integral, passing to gauge pressure
\begin{align}
\int_{P_s(o,t)}^{P_s(L,t)}(P_s(x,t)+b) \ dP_s(x,t).
\end{align}
Since $P_s=P_g+P_a$, all you need to do is that change of variables. Besides, $dP_s=dP_g$ because $P_a$ is a constant. Also, there is no need to say $dP_s(x,t)$, since these are dummy variables.
\begin{align}
\int_{P_s(o,t)}^{P_s(L,t)}(P_s+b) \ dP_s&=\int_{P_s(o,t)}^{P_s(L,t)}(P_g+P_a+b)dP_s\\
&=\int_{P_s(o,t)}^{P_s(L,t)}(P_g+P_a)dP_s+(P_{s,L}-P_{s,o}) b\\
&=\int_{P_s(o,t)}^{P_s(L,t)}P_g dP_s+(P_{s,L}-P_{s,o}) (b+P_a)\\
&=\int_{P_g(o,t)}^{0}P_g dP_g-P_g (b+P_a)\\
&=-\frac{P_g^2}{2}-P_g (b+P_a)=-\frac{P_g}{2}(P_g+2P_a+2b),\\
\end{align}
where in the last line $P_g$ is $P_{g,o}(t)$. Note that $P_{s,L}-P_{s,o}=-P_g$.
How does dividing atmospheres by absolute pressure give you gauge pressure?
Here I think you are confused. psi, psig, atm, and Pa are all units of pressure. That is why you can add gauge pressure to atmospheric pressure in the first place. The problem is that pressures expressed in these units are not proportional since they have different zero, which is the same problem when working between the different temperature scales (farenheit, celsius, and Kelvin). My advice is to make every algebraic and mathematical manipulation without units, or using a single system like SI, and only translate to more practical units at the end.
If you neglect viscosity, Bernoulli's equation (just Navier-Stokes without frictional or stress terms) will get you into the ballpark:
$$P_g + \frac{1}{2}\rho_g v_g^2 = P_a$$
Where the $g$ subscripts pertain to the gas and the $a$ subscript to the ambient. The gas density $\rho_g \equiv M / V$ is the ratio of the mass of gas (M) in the tank to the volume of the tank. If the tank is a rigid container (like a propane tank) then the volume of the gas is constant and the pressure will vary with the mass flow and the temperature. If you assume the tank remains at a constant ambient temperature, the pressure will only vary with the mass flow rate (isothermal expansion) and you can obtain that from the ideal gas law:
$$P_g = \frac{m}{M} RT$$
where $m$ is the molecular mass of the gas in question, $T$ the temperature, $R$ the gas constant, and $M$ the total mass of the gas remaining in the tank. This is a function of time because mass is leaving the tank. The rate at which mass leaves is a function of the exit velocity (it depends on the volumetric flow rate, which is a product of the exit orifice size and the exit velocity). So you can solve for $M(v_g)$ and substitute in the above equations and solve for $v_g$ self-consistently. Note this approach also ignores any pipes that might be attached to the orifice. For that, you'd need to calculate the volumetric flow rate using Poiseuille's equation.
Best Answer
Starting with the conservation of mass for the system
$$ \frac{dm}{dt} = -\rho v A \tag{1}$$
with $m$ the mass of the gas, $\rho$ its density, $v$ its velocity and $A$ the area of the hole
where $m = n/M_w$ and $ n = PV/RT $, this simplifies to
$$ \frac{d}{dt}(\frac{PV}{RTM_w}) = -\rho v A \tag{2}$$
The gas volume will remain constant, but the pressure will not (im assuming constant T as well, it gets real messy if its not ). Simplifying a bit substituting the density with the ideal gas law, we find
$$ \frac{V}{RTM_w}\frac{dP}{dt} = -\frac{M_wP}{RT} v A \tag{3}$$
Simplifying a bit
$$ \frac{1}{P}\frac{dP}{dt} = -\frac{M_w^2A}{V} v \tag{4}$$
You wish to find $P$ but cant integrate just yet because $v$ is time and pressure dependent, so we need expressions for $v$
So, according to Bernoulli,
$$v^2 = \frac{2(P-P_{atm})}{\rho} \tag{5}$$
again, replacing the density with ideal gas law,
$$v^2 = \frac{2(P-P_{atm})RT}{M_wP} \tag{6}$$
Simplifying a bit
$$v^2 = \frac{2RT}{M_w} - \frac{2RTP_{atm}}{M_w}\frac{1}{P} \tag{7}$$
Taking the square root and grouping constants,
$$v = \sqrt{\alpha - \beta\frac{1}{P}} \tag{8}$$
Substituting (8) into (4) and grouping the constants in (4), we end up with
$$\frac{dP}{dt} = -\gamma P \sqrt{\alpha - \beta\frac{1}{P}} \tag{9}$$
Now you can integrate. That is waaaay above my skill level to do analytically but wolfram did something
$$P(t) = \frac{e^{-c_1 \sqrt{\alpha} -\gamma t \sqrt{a}} (\beta e^{\gamma t\sqrt{\alpha}} + e^{c_1 \sqrt{a} })^2}{4 a}$$
No idea if that is even remotely correct. I would just do equation (9) numerically
EDIT: Did it numerically with your data
I cant figure how to change the colour of the graph in Polymath but I think its good enough to get a general idea. The gas losses 99.99698% of its initial pressure over 150 seconds. The wolfram analytical answer might actually be correct it seems
EDIT2: In your practical experiments it might take longer because of cooling effects of the gas, but I dont know how much longer. A good starting point if you want to derive such a system, is to start at equation (2) and use the chain rule on the P/T terms (Volume of the system is constant, so isochoric equations might help). It produces the problem of having to find equations for P' and v as before, but also T'. And I dont have a plan to do that at the moment because I cant seem to use the first law of thermo without producing another term (Q' in this case) because the PV work is zero, and therefor have no way of relating T' to anything else.
EDIT3: Perhaps one could use an enthalpy approach. $dH = TdS + VdP$ and $dS = dQ/T$ , so that $dH = dQ + VdP$.
Therefor $C_pdT = C_vdT + VdP$ Simplifying delivers $(C_p-C_v)dT = VdP$
Here is where it starts getting tricky again. Normally, Cp and Cv are temperature dependent. Divide them out, and then you have an equation for dT (albeit dirty) that can be replaced in the modified equation (2) so that equation (2) is only an ODE in Pressure. Solving that one and the temperature one simultaneously might produce a more realistic answer. I could be wrong though, so take EDIT3 with a grain of salt
EDIT4: Reply to comments got too long
Changing it to allow for a slightly larger initial pressure is a trivial exercise, really. Graph shape will remain the same, initial condition only changes. Regarding the "realism" of it, it might be faster in real life because a quick check on the gas velocity showed that it was just over the speed of sound (initially). This means your well above the rule of thumb of mach 0.3 limit for in-compressible fluids. I do not have enough practical or mathematical experience to give an estimate on how a real gas would change the system. Doing that is a lot of work, where your choice of assumptions would definitely impact the result, e.g fugacity based approach. For your second comment, yes of course. Simply use the ideal gas law to relate the quantities. $\rho = PM_w/RT$. Doing it in chunks, yes of course. I manually chose the run time of the model. In chunks you will need better software than Polymath so that you have exit conditions and such for your loops. Matlab would be great at it. Simulink would work too if you know/can find out how triggers work. If you want to iterate every 10psi to find a "simpler" solution to the problem, you will merely have an approximate version of what I have here. Here, essentially, Im running it in 0.002psi increments.