Let $G$ be a locally compact Hausdorff group with a left Haar measure $\lambda$. Define the convolution of two functions $f,g \in L^1(G)$ by
$$(f \ast g)(x) = \int f(y) g(y^{-1}x) d\lambda (y), ~~~ \forall x \in G$$
If the group $G$ is abelian the convolution is commutative: $f \ast g = g \ast f$.
In general, for any $x \in G$ we have (written multiplicatively)
$$ (f \ast g)(x) = \int f(y) g(y^{-1}x) d\lambda(y) = \int f(xy) g((xy)^{-1}x) d\lambda(y) = \int f(xy) g(y^{-1}) d\lambda(y)$$
In the second equality, we apply a left shift by $x^{-1}$ which does not change the integral since $\lambda$ is left invariant.
Precomposing with inversion yields
$$ \int f(xy^{-1}) g(y) d\rho(y)$$
where $\rho$ is the associated right Haar measure defined by $\rho(B) = \lambda(B^{-1})$ for any Borel set $B \subseteq G$.
Finally, commuting $x$ and $y^{-1}$ gives
$$ \int g(y) f(y^{-1}x) d\rho(y)$$
Now, if $G$ is unimodular, $\rho$ and $\lambda$ coincide, so the last expression is the convolution $g \ast f$. Also, since both $y^{-1} \in G$ and $x \in G$ are arbitrary, the step requires $G$ to be abelian (which then also makes it unimodular).
I am looking for an explicit counterexample to the claim that $f \ast g = g \ast f$ in general, and conditions under which the formula is true (which are hopefully weaker than $G$ being abelian).
Thank you very much in advance!
As you noted if $G$ is abelian then it is trivial that convolutions commute.
For the converse, let convolution of any two $C_c$ functions commute. Let $f,g \in C_c(G)$
Then $\forall x \in G \text{ we have }$ $$0= f*g(x)-g*f(x)=\int_G f(xy)g(y^{-1}) d\lambda(y) - \int_{G} g(y)f(y^{-1}x)d\lambda(y)$$ $$=\int_G f(xy^{-1})g(y)\Delta(y^{-1}) d\lambda(y) - \int_{G} g(y)f(y^{-1}x)d\lambda(y)$$ $$\implies \int_G g(y)(\Delta(y^{-1})f(xy^{-1})-f(y^{-1}x))d\lambda(y)=0$$
Since, $g \in C_c(G)$ was arbitrarily chosen, it follows that $$\Delta(y^{-1})f(xy^{-1})=f(y^{-1}x), \forall x,y \in G$$ So put $x=1$ above and note that $\Delta(y^{-1})f(y^{-1})=f(y^{-1})$ . Again $f \in C_c(G)$ was arbitrarily chosen thus $f$ can very well be non-zero at $y^{-1}$. So we get, $\Delta(y^{-1})=1, \forall y \in G$
Hence, $f(xy^{-1})=f(y^{-1}x) \forall x,y \in G$ . Then just replace $y$ by $y^{-1}$ and we get $$f(xy)=f(yx) \forall f \in C_c(G) \implies xy=yx, \forall x,y \in G$$
Since, you have the result for $C_c(G)$, it follows for $L^1(G)$