When dealing with any mathematical problem, often the best way to do it is to break it up into smaller parts that are easier to deal with. So starting with our original linear inhomogeneous IVP,
$$(\partial_t-\partial_x^2)u(t,x)=h(t,x) \\ u(0,x)=u_0(x) \\ (t,x)\in\mathbb R_{\geq 0}\times \mathbb R \tag{0}$$
We consider the simpler problems
$$(\partial_t-\partial_x^2)u(t,x)=0 \\ u(0,x)=u_0(x)\tag{1}$$
$$(\partial_t-\partial_x^2)u(t,x)=h(t,x) \\ u(0,x)=0 \tag{2}$$
(And, in order for uniqueness, all problems are coupled with the additional physical constraint that the integral of $|u|^2$ along the $x$ axis is finite for all time)
If we have found a solution $u_1$ to (1) and a solution $u_2$ to (2) then due to the linearity of the heat operator, $u_1+u_2$ is a solution of (0).
What is nice about this approach is that (1) and (2) are already well studied problems. For the first, we use the fundamental solution of the heat equation, also known as the heat kernel, and for the second, we use the Green's function (see item 14 in the linked table).
In the following, assume $(\tau,\xi)\in\mathbb{R}_{\ge 0}\times\mathbb R$.
The fundamental solution is the solution $F(t,x;\xi)$ of the homogeneous IVP
$$(\partial_t-\partial_x^2)u(t,x)=0 \\ u(0,x)=\delta(x-\xi)$$
This is a very well known problem, and the solution is
$$F(t,x;\xi)=\frac{\exp\left(-\frac{(x-\xi)^2}{4t}\right)}{\sqrt{4\pi t}}$$
The Green's function is the solution $G(t,x;\tau,\xi)$ of the inhomogeneous IVP with zero initial conditions,
$$(\partial_t-\partial_x^2)u(t,x)=\delta(t-\tau)\delta(x-\xi) \\ u(0,x)=0$$
And again, the solution of this problem is also well known,
$$G(t,x;\tau,\xi)=\frac{\exp\left(-\frac{(x-\xi)^2}{4(t-\tau)}\right)}{\sqrt{4\pi (t-\tau)}}=F(t-\tau,x;\xi)$$
Remark: The curious, and very non-obvious symmetry between the fundamental solution and the Green's function is an example of Duhamel's Principle, and applies for a broad range of linear evolution equations, i.e equations of the form
$$\partial_t u=\mathscr L [u]$$
Where $\mathscr L$ is some linear differential operator involving only $u$ and its spacial derivatives.
Now that these two solutions are known, we can write $u_1$ and $u_2$ as convolutions,
$$u_1=F\ast_x u_0 \\ u_2 = G\ast_{(t,x)} h$$
I.e,
$$u_1(t,x)=\int_{\mathbb{R}} F(t,x;\xi) u_0(\xi)\mathrm d\xi \\ u_2(t,x)=\int_0^t \int_{-\infty}^\infty G(t,x;\tau,\xi)h(\tau,\xi)\mathrm d\xi \mathrm d\tau$$
And,
$$u=u_1+u_2$$
In your case with $u_0(x)=\cos(3x)~,~h(t,x)=e^t$, you end up with
$$u(t,x)=\int_{\mathbb{R}} \frac{\exp\left(-\frac{(x-\xi)^2}{4t}\right)}{\sqrt{4\pi t}} \cos(3\xi)\mathrm d\xi+\int_0^t \int_{-\infty}^\infty \frac{\exp\left(-\frac{(x-\xi)^2}{4(t-\tau)}\right)}{\sqrt{4\pi (t-\tau)}}e^\tau \mathrm d\xi \mathrm d\tau$$
Though often this procedure will not give you closed forms all that easily, it is the best way to deal with these kinds of problems in general, because typically closed forms are just not possible. To get closed forms, one would usually use separation of variables, but I stress that, in general, techniques like separation of variables will just fail.
For the relevant theory on the above work, see chapter 9 of Peter J. Olver's Introduction to Partial Differential equations, in particular pages 292 to 297.