We see that $$f(t) = e^{-3t} + e^{-t} \int_0^t e^{-\tau}f(\tau) d\tau$$
I am having issues solving for $f(t)$. I recognize that if we set $g(\tau) = e^{-\tau} f(\tau)$ and if we set $h(t) = \int_0^t g(\tau) d\tau,$ Then $$f(t) = e^{-3t} + e^{-t}h(t)$$ $$\implies \mathscr{L}\{f(t)\} = \mathscr{L}\{e^{-3t} + e^{-t}h(t)\} = \mathscr{L}\{e^{-3t}\} + \mathscr{L}\{e^{-t}h(t)\}$$ $$= \frac{1}{s+3} + H(s + 1) = \frac{1}{s+3} + \mathscr{L}\{h(t)\}\mid_{s \rightarrow s+1}$$ $$= \frac{1}{s+3} +\mathscr{L}\{\int_0^t g(\tau) d\tau\}\mid_{s \rightarrow s+1} = \frac{1}{s+3} + \frac{G(s + 1)}{s+1}.$$
However, I feel like I'm stuck here, because when I take the inverse Laplace Transform I end up this what I started with. Did I miss a step? I feel like I should be doing something with the convolution theorem.
Let $g(s)=(\mathcal{L} f)(s)$. The Laplace transform of $e^{-x}f(x)$ is given by $g(s+1)$ and the Laplace transform of $\int_{0}^{x}e^{-u}f(u)\,du$ is given by $\frac{1}{s}g(s+1)$, hence your differential equation can be written in terms of $g$ as
$$ g(s) = \frac{1}{s+3}+\frac{g(s+2)}{s+1}\tag{1} $$ leading to: $$ g(s) = \frac{1}{s+3}+\frac{1}{s+1}\left(\frac{1}{s+5}+\frac{1}{s+3}\left(\frac{1}{s+7}+\frac{1}{s+5}\left(\frac{1}{s+9}+\ldots\right)\right)\right)\tag{2}$$
$$\small g(s)=\frac{1}{s+3}+\frac{1}{(s+1)(s+5)}+\frac{1}{(s+1)(s+3)(s+7)}+\frac{1}{(s+1)(s+3)(s+5)(s+9)}+\ldots $$ Now I suggest the following approach: $g(s)$ is a meromorphic function with poles at the negative odd integers. It is simple to compute the closed form of $R_\xi=\text{Res}_{s=\xi}\,g(s)$ for any $\xi\in\Xi=\{-1,-3,-5,\ldots\}$, then to consider the inverse Laplace transform of $g$:
$$ f(x) = \sum_{\xi\in\Xi}R_\xi e^{\xi x}.\tag{3} $$ Mathematica finds a nice closed form for $g(s)$: $$\begin{eqnarray*} g(s) &=& \frac{1+2^{\frac{1+s}{2}} \sqrt{e} \left(-\text{Gamma}\left[\frac{3+s}{2}\right]+\text{Gamma}\left[\frac{3+s}{2},\frac{1}{2}\right]\right)}{1+s}\\&=&\frac{1}{1+s}\left[1-2^{\frac{1+s}{2}}\sqrt{e}\int_{0}^{1/2}u^{\frac{s+1}{2}}e^{-u}\,du\right]\\&=&\frac{2-\sqrt{e}}{s+1}+\frac{\frac{1}{2}\sqrt{e}}{s+3}-\frac{\frac{1}{8}\sqrt{e}}{s+5}+\frac{\frac{1}{48}\sqrt{e}}{s+7}-\frac{\frac{1}{384}\sqrt{e}}{s+9}+\ldots\tag{4} \end{eqnarray*}$$ and from the last line it is straightforward to recover a solution $f(x)$:
$$\begin{eqnarray*} f(x)&=&(2-\sqrt{e})e^{-x}+\sqrt{e}\sum_{n\geq 1}\frac{(-1)^{n+1}}{(2n)!!} e^{-(2n+1)x}\\&=&\color{red}{e^{-x} \left(2-e^{e^{-x} \sinh(x)}\right)}.\tag{5}\end{eqnarray*} $$