The complete list of irreducible admissible representations of $\mathrm{GL}(2, \mathbb{R})$ are known - principal series, discrete series, limit of discrete series, and finite dimensional representations (symmetric power of standard representation). In Bump's book, chapter 2.8, he introduces the Whittaker model and shows that if the Whittaker model exists, then unique.
However, it seems that he doesn't talk about existence at all. (He talked about the existence of Whittaker models for infinite dimensional representations over non-archimedean local fields later, by using the theory of Jacquet functor).
I believe that most of the irreducible admissible representations have Whittaker models, but I have no idea how to construct it.
For example, assume that we have a principal series $\mathcal{P}_{\mu}(\lambda, \epsilon)$, where $\mu, \lambda\in \mathbb{C}$ are scalars corresponds to the elements $Z = \begin{pmatrix}1&0\\0&1\end{pmatrix}, \Delta\in Z(U\mathfrak{gl}(2, \mathbb{R}))$, and $\epsilon\in \{0, 1\}$ represents its parity. Assume that it is irreducible (i.e. $\lambda \neq \frac{k}{2}(1-\frac{k}{2})$ for all $k\equiv \epsilon$ (mod 2). How can we construct Whittaker model for this?
For any local field $k$, an intertwining operator from the induced model for a principal series to the Whittaker model (such as referred to by Paul Garrett in his comment above) is explicitly: $$\phi \mapsto W_\phi,$$ where $$W_\phi(g) = \int_{k} \phi \left( ( \begin{smallmatrix} & -1 \\ 1 & \end{smallmatrix}) ( \begin{smallmatrix} 1 & x\\ & 1 \end{smallmatrix}) g \right) \psi(x)\,dx,$$ where $\psi$ is the chosen additive character for your Whittaker model. In general, this only converges conditionally, but if you interpret the integral as the limit over an increasing sequence of compact subsets of $k$ defined by $|x|\leq R$, $R \to \infty$, then it does give an intertwiner.