How to solve $\nabla^2\phi(\textbf{r})=m^2\phi(\textbf{r})$? How can I solve an equation of the form $\nabla^2\phi(\textbf{r})=m^2\phi(\textbf{r})$ where $m^2$ is a real positive constant and $\textbf{r}=(x,y,z)$ is a point 3-D space. I want to write the solution at least formally and Fourier transform doesn't help. The Fourier transform $$\phi(\textbf{r})=\int\frac{d^3\textbf{k}}{(2\pi)^{3/2}}e^{i\textbf{k}\cdot\textbf{r}}\tilde{\phi}(\textbf{k})\tag{1}$$ gives $$(\textbf{k}^2+m^2)\tilde{\phi}(\textbf{k})=0\tag{2}.$$I want to find a nonzero solution $\phi(\textbf{r})$ from (1) with the constraint $k^2+m^2\neq 0$. I have no idea how to proceed next.

The only physical solution here is $\phi \equiv 0$ as your equation says that $\phi({\bf k})=0$ since $k^2+m^2\not = 0$. Note that this derivation assumes that the Fourier transform of $\phi$ exists (which is in most physical applications is a very reasonable assumption).
Another way to derive this: if we have a general region $V$ and impose the boundary condition $\phi = 0$ on $\partial V$ then by multiplying the PDE by $\phi$ and integrating using the we get 
$$0 = \int_V m^2 \phi^2 - \phi\nabla^2\phi\,{\rm d}{\bf x} = \int_V (\nabla \phi)^2 + m^2 \phi^2\,{\rm d}{\bf x}$$
and both terms are strictly positive which means that $\phi \equiv 0$ in $V$. In your case $V = \mathbb{R}^3$ and $\partial V$ is "at infinity" and I'm assuming here $\phi$ decays "fast enough" to $0$ as ${\bf r}\to \infty$ as for the integral above to exist. The integral on the right hand side above is what we usually call the "energy" of the field so if any other solution exists then it must correspond to infinite "energy".
There are infinitely many solutions of this kind: $\phi({\bf r}) = c\frac{\sinh(mr)}{r}$ where $r = \|{\bf r}\|$ and $\phi({\bf r}) = Ae^{ax + by + cz}$ if $a^2+b^2+c^2 = m^2$ are two simple examples. What breaks down in our original derivation above for these solutions is that the Fourier transform simply does not exist.
$$ $$
