We are given that \( f(x + y) = f(x) f(y) \), which implies that \( f(x) = e^{\lambda x} \), for some constant \( \lambda \). We also know that \( f'(0) = 4a \), so:
\[
f'(x) = \lambda e^{\lambda x}, \quad f'(0) = \lambda = 4a.
\]
Thus, \( f(x) = e^{4ax} \).
Next, we are given the differential equation \( f''(x) - 3a f'(x) - f(x) = 0 \), and we need to solve for \( a \). Substituting \( f(x) = e^{4ax} \) into the equation gives:
\[
f''(x) = 16a^2 e^{4ax}, \quad f'(x) = 4a e^{4ax}.
\]
Substitute into the differential equation:
\[
16a^2 e^{4ax} - 3a (4a e^{4ax}) - e^{4ax} = 0,
\]
\[
16a^2 - 12a^2 - 1 = 0 \quad \Rightarrow \quad 4a^2 = 1 \quad \Rightarrow \quad a = \frac{1}{2}.
\]
Thus, the function becomes:
\[
f(x) = e^{2x}.
\]
The area of the region is given by:
\[
\text{Area} = \int_0^2 f(ax) \, dx = \int_0^2 e^{x} \, dx = e^2 - 1.
\]
Thus, the area of the region is \( \boxed{e^2 - 1} \).