Hmm, after a bit of googling I found that my conjecture is in fact proved in the 2018 Harvard honors thesis "Modular Magic" by Aaron Slipper. (The lucky google search that yielded this was "2607840 theta series" -- 2607840 being the coefficient of $q^6$ in $f(q)/2$ -- but it only occurred to me to try this after an evening spent googling and reading related literature.)
Funnily enough, Slipper proved the fact that I need as part of his exposition of Viazovska's work on sphere packing, and needed it for reasons that seem a bit related to my own reasons for needing it, but not exactly to prove Viazovska's inequalities in the simple fashion I discovered.
Slipper's argument, found in page 76 of his thesis (as a component of his proof of Proposition 4.4.6 from page 75), is quite elegant. In my notation, he found that the function $f$ can be represented as
$$
f(q) = \frac{2 \Theta(q) + \Theta(-q)}{3}
$$
where $\Theta(q)$ is the theta function associated with a certain 20-dimensional lattice, called the DualExtremal(20,2)a lattice. Since the Fourier coefficients of the theta series of a lattice are nonnegative by definition, we immediately get the claim about the nonnegativity of the coefficients of $f(q)$.
Now, this is already very nice, but I'm wondering if this argument can be reworked to get a direct proof of the nonnegativity from simple properties of the Jacobi thetanull functions. In particular, I worked out that the function $\Theta(q)$ can be represented as
$$
\Theta = 2X^5 - 5 X^4 Y + 5 X^3 Y^2 + 5 X^2 Y^3 - 5 X Y^4 + 2Y^5,
$$
where $X=\theta_3^4$, $Y=\theta_2^4$. So maybe from here it's not too hard to show that $\Theta$ has nonnegative Fourier coefficients. I'll think more about this tomorrow. Edit: yes, I was right about that -- see below.
Edit: here is a self-contained proof of my nonnegativity conjecture, inspired by the argument from Slipper's thesis but reworked to give a more direct proof that does not rely on facts about lattices and their theta series:
Denote $X=\theta_3^4$, $Y=\theta_2^4$ as above.
Because of the well-known identity $\theta_3^4=\theta_2^4+\theta_4^4$, $f(q)$ can be written as
\begin{align*}
f &= X^3 Y^2 + X^2 Y^3 + X^3 (X-Y)^2 + X^2 (X-Y)^3
\\ &= 2 X^5 - 5 X^4 Y + 5 X^3 Y^2
\end{align*}
Define a function $\Theta(q)$ by
$$
\Theta(q) = 2X^5 - 5 X^4 Y + 5 X^3 Y^2 + 5 X^2 Y^3 - 5 X Y^4 + 2Y^5.
$$
Observe that, since $X(-q)=X-Y$ and $Y(-q)=-Y$, we have
\begin{align*}
\Theta(-q) &= 2 (X-Y)^5 - 5 (X-Y)^4 (-Y) + 5 (X-Y)^3 (-Y)^2
\\ & \qquad + 5 (X-Y)^2 (-Y)^3 - 5 (X-Y) (-Y)^4 + 2 (-Y)^5
\\ &=
2 X^5 - 5 X^4 Y + 5 X^3 Y^2 - 10 X^2 Y^3 + 10 X Y^4 - 4 Y^5.
\end{align*}
From the last two equations it follows that
$$
f(q) = \frac{2\Theta(q) + \Theta(-q)}{3}.
$$
Note that an alternative way to express $\Theta$ is as
$$
\Theta = \frac{(2X-Y)^5 + 30 (2X-Y)^2 Y^3 + 15 (2X-Y) Y^4 + 18 (2X-Y)^5}{16}
$$
Note that both $Y$ and $2X-Y$ have nonnegative Fourier coefficients. For $Y$ this is immediate from the definition, and for $2X-Y$ this follows from the identity
$$
\frac{2X-Y}{2} = 1 + 24 \sum_{n=1}^\infty \sigma_{\textrm{odd}}(n) q^{2n},
\qquad (*)
$$
where $\sigma_{\textrm{odd}}(n)$ is the odd divisor function, defined by
$$
\sigma_{\textrm{odd}}(n) = \sum_{d\,|\,n,\ \ d\textrm{ odd}} d.
$$
We expressed $\Theta$ as a linear combination with positive coefficients of monomials in $Y$ and $2X-Y$, which have nonnegative Fourier coefficients. Therefore $\Theta$ has nonnegative Fourier coefficients.
By the way we expressed $f$ in terms of $\Theta$, it's clear that $f$ too must have nonnegative Fourier coefficients. QED.
Remarks:
The identity $(*)$ is standard, see for example this OEIS page. (I haven't checked, but I think it can be proved using the identities in this MathOverflow post expressing $\theta_2^4$ and $\theta_3^4$ in terms of Eisenstein series.)
The modular form $2X-Y$ is twice the Eisenstein series of weight 2 referred to in Henri Cohen's answer, where it was written as $T^4-8F_2$. (Thanks Henri for the suggestion to look for a representation of $f$ involving this function!)
The function I'm denoting by $\Theta$ is the same as the function $\Theta$ from page 76 of Slipper's thesis. But I defined it in terms of the Jacobi thetanulls $\theta_2$, $\theta_3$, and he defined his $\Theta$ as the theta series of a certain (pretty esoteric as far as I can tell) 20-dimensional lattice. Proving that the two functions coincide would require additional work, but isn't necessary to make my version of the nonnegativity proof work.
The proof above does not explain how I found the representation of $\Theta$ as a linear combination with positive coefficients of monomials in $2X-Y$ and $Y$. I did this by solving a linear program in 21 variables (using Mathematica's Minimize
command), but maybe there's an easier way.