Looking at the table here you will recognize three different possible behaviors. Let us see why. Consider the denominator. This is can be rewritten as
$$s^2+\frac{R}{L}s+\frac{1}{LC}=(s+\alpha)^2+\beta^2$$
where
$$\alpha=\frac{R}{2L} \qquad \beta=\sqrt{\frac{1}{LC}-\frac{R}{2L}}.$$
So, when $\frac{1}{LC}>\frac{R}{2L}$ you will recognize an exponentially decaying sine wave. When $\frac{1}{LC}=\frac{R}{2L}$ you will get just an exponential decay. When $\frac{1}{LC}<\frac{R}{2L}$ you will get an exponential decay multiplied by a hyperbolic cosine. All this can be deduced from the table I linked at the beginning of this answer.
I realized that the easier approach than trying to evaluate what I had reached in my derivation is to introduce some dimensionless variables to simplify things.
We have just like above that
$$
\frac{dm}{dt} = q\cdot \alpha\left(1 - H(t-T)\right) - q \cdot \frac{m}{V_{0}}.
$$
Then by dimensional analysis, we have that $\left[\frac{V_{0}}{q}\right] = T$, and so we introduce a characteristic time by $t_{c} = \frac{V_{0}}{q}$. Then $\tau = \frac{t}{t_{c}}$ is dimensionless and by the Chain rule $\frac{dm}{dt} = \frac{dm}{d\tau}\frac{d\tau}{dt} = \frac{q}{V_{0}}\frac{dm}{d\tau}$. We also introduce a dimensionless mass $\mathcal{M}$ via $m_{c} = V_{0}\alpha$, and upon subbing dimensionless quantities into the DE have
$$
\frac{d\mathcal{M}}{d\tau} = 1 - H\left(\tau - \frac{qT}{V_{0}}\right).
$$
From here, the Laplace transform is easily performed to acquire
$$
\mathcal{L}[\mathcal{M}] = \frac{1}{s} - \frac{1}{s+1} - e^{\frac{-qTs}{V_{0}}}\left(\frac{1}{s} - \frac{1}{s+1}\right).
$$
This is much easier to apply the inverse Laplace operator to than what I had in my original question, and using the shift theorems
$$
\mathcal{L}^{-1}\left[\mathcal{L}\left[\mathcal{M}\right]\right] = \mathcal{M} = 1 - e^{-\tau} - H\left(\tau - \frac{qT}{V_{0}}\right)\left[1 - e^{-\left(\tau - \frac{qT}{V_{0}}\right)}\right].
$$
From here, substituting dimensional values back in and evaluating at $t = 2T$ isn't too hard.
Best Answer
When all of the coefficients are either linear functions or constant functions, there is a so call "kernel method" which is modified from the Laplace transform method to solve them:
Let $y=\int_Ce^{xs}K(s)~ds$ , then you can obtain the first-order ODE for $K(s)$ . Once you find $K(s)$ , you can substitute back to $\int_Ce^{xs}K(s)~ds$ . Since the procedure in fact suitable for any complex number $s$ , you can rewrite the integral form as $\int_{a_n}^{b_n}e^{x(p_n+q_ni)t}K((p_n+q_ni)t)~d((p_n+q_ni)t)$ . By choosing suitable $x$-independent real numbers groups of $a_n$ , $b_n$ , $p_n$ and $q_n$ , you can find the linear independent solutions in integral forms.
For example in How do you solve $y''+\frac{x}{2}y'+y=0$ a 2nd order homogenous equation?.
Sometimes you need to find the linear independent solutions with the combinations of different types of integral forms , for example in Particular solution to a Riccati equation $y' = 1 + 2y + xy^2$.