How to solve this Poisson's equation I want to solve Poisson's equation given by
$$
\partial_t f = \nabla^2f = \partial^2_r f + \frac{1}{r}\partial_r f
$$
on a circular disc, so it is 2D. I use separation of variables to write $f=R(r)T(t)$, which gives me these equations for $R$ and $T$
$$
T'+\lambda T=0\\
rR'' + R' + \lambda r R=0
$$
The equation for $T$ is trivial to solve, but how do I solve the one for $R$?
 A: Your equation is:
$$r{\frac {d^{2}}{d{r}^{2}}}R \left( r \right) +{\frac {d}{dr}}R \left( 
r \right) +\lambda\,rR \left( r \right)=0 \tag{1}$$ 
if we multiply $(1)$ by $r$ and introduce the function:
$$S(\sqrt{\lambda}\,r)=R(r)\tag{2}$$
then $(1)$ becomes:
$$r^2{\frac {d^{2}}{d{r}^{2}}}S(\sqrt{\lambda}\,r) +r{\frac {d}{dr}}S(\sqrt{\lambda}\,r) +\lambda\,r^2S(\sqrt{\lambda}\,r)=0 \tag{3}$$ 
if we then change variable by making the substitution $$s=r\,\sqrt{\lambda} \tag{4}$$
we obtain (see Appendix A):
$$s^2{\frac {d^{2}}{d{s}^{2}}}S \left( s \right) +s{\frac {d}{ds}}S \left( 
s \right) +s^2S \left( s \right)=0 \tag{5}$$ 
which is a Bessel differential equation and its solution is:
$$S \left( s \right) =A\,
{{\rm J}_0\left(s\right)}+B\,
{{\rm Y}_0\left(s\right)} \tag{6}$$
hence:

$$R \left( r \right) =S(\sqrt{\lambda}\,r)=A\,
{{\rm J}_0\left(\sqrt {\lambda}r\right)}+B\,
{{\rm Y}_0\left(\sqrt {\lambda}r\right)} \tag{7}$$

Constants $A$ and $B$ would then be fixed subject to boundary conditions, for example, $B$ may be zero as ${{\rm Y}_0\left(\sqrt {\lambda}r\right)}$ has a singularity at $r=0$. 
Appendix A
To convince ourselves of the variable change from $(3)-(5)$ consider the limit definition of the derivative. After multiplying $(1)$ by $r$ we have:
$${r}^{2}\lim _{h\rightarrow 0}{\frac {R \left( r+h \right) -2\,R
 \left( r \right) +R \left( r-h \right) }{{h}^{2}}}+r\lim _{h
\rightarrow 0}{\frac {R \left( r+h \right) -R \left( r \right) }{h}}+
\lambda\,{r}^{2}R \left( r \right) =0$$
and according to $(2)$ we have:
$$R \left( r \right) =S \left( \sqrt {\lambda}\,r \right) ,\quad R
 \left( r-h \right) =S \left( \sqrt {\lambda}\,r-\sqrt {\lambda}\,h
 \right) ,\quad R \left( r+h \right) =S \left( \sqrt {\lambda}\,r+\sqrt {
\lambda}\,h \right)  $$
hence:
$${r}^{2}\lim _{h\rightarrow 0}{\frac {S \left( \sqrt {\lambda}r+\sqrt {
\lambda}h \right) -2\,S \left( \sqrt {\lambda}r \right) +S \left( 
\sqrt {\lambda}r-\sqrt {\lambda}h \right) }{{h}^{2}}}+r\lim _{h
\rightarrow 0}{\frac {S \left( \sqrt {\lambda}r+\sqrt {\lambda}h
 \right) -S \left( \sqrt {\lambda}r \right) }{h}}+\lambda\,{r}^{2}S
 \left( \sqrt {\lambda}r \right) =0$$
after the variable change $(4)$ this becomes:
$${s}^{2}\lim _{h\rightarrow 0}{\frac {S \left( s+\sqrt {\lambda}h
 \right) -2\,S \left( s \right) +S \left( s-\sqrt {\lambda}h \right) }
{\lambda\,{h}^{2}}}+s\lim _{h\rightarrow 0}{\frac {S \left( s+
\sqrt {\lambda}h \right) -S \left( s \right) }{\sqrt {
\lambda}\,h}}+{s}^{2}S \left( s \right) =0$$
after introducing the new infinitesimal $\epsilon=\sqrt{\lambda}\,h$ this becomes:
$${s}^{2}\lim _{\epsilon\rightarrow 0}{\frac {S \left( s+\epsilon
 \right) -2\,S \left( s \right) +S \left( s-\epsilon \right) }{{
\epsilon}^{2}}}+s\lim _{\epsilon\rightarrow 0}{\frac {S \left( s+
\epsilon \right) -S \left( s \right) }{\epsilon}}+{s}^{2}S \left( s
 \right) =0$$
and after taking the limit $(5)$ follows.
Appendix B
Update:
To compare with @user8268's method note that:
$$z=\sqrt{\lambda}\,y:\quad\int _{-\sqrt {-\lambda}}^{\sqrt {-\lambda}}\!{\frac {{{\rm e}^{rz}}}{
\sqrt {{z}^{2}+\lambda}}}{dz}=i\int _{-1}^{1}\!{\frac {{{\rm e}^{ir
\sqrt {\lambda}y}}}{\sqrt {1-{y}^{2}}}}{dy}=i\pi \,{{\rm J}_0\left(r\sqrt {\lambda}\right)}$$
by the Poisson representation of the Bessel function.
A: The equation for $R$ is a "linear differential equation with linear coefficients", i.e. of the form
$$\sum_{k=1}^n(a_k\,r+b_k)\frac{d^k}{dr^k}R=0\qquad(*).$$
There is a quite beautiful way for solving these equations, especially for geting solutions with nice asymptotic properties. I'll describe it, though it's a bit overkill just to solve your single equation. Several interesting equations (e.g. the Schrodinger equation for the hydrogen atom) reduce to the form $(*)$, so it's quite useful to know this method.
The idea is roughly that Fourier (or Laplace) transform exchanges derivatives with multiplication by the variable, so it will change $(*)$ to a 1st order equation.
Let $A(z)=\sum_k a_k z^k$ and $B(z)=\sum_k b_k z^k$. Let us look for a solution of $(*)$ of the form
$$R(r)=\int_C F(z)e^{rz}dz\qquad(\times)$$
where $F$ is a function of complex variable $z$ and $C$ is an appropriate path in $\mathbb C$. The equation $(*)$ becomes (if we can exchange integral and derivatives)
$$\int_C(A(z)r+B(z))F(z)e^{rz}dz=0.$$
This is satisfied in the case when 
$$(A(z)r+B(z))F(z)e^{rz}=d(G(z)e^{rz})/dz\qquad(+)$$
 where $G$ is another holomorphic function, provided that $G(z)e^{rz}$ goes to (say) $0$ at the endpoints of $C$ (for $r$'s on the interval where you want to solve $(*)$, in your case $(0,\infty)$), or if $C$ is closed (and $G$ is univalued along $C$).
Equation $(+)$ is $A(z)F(z)=G(z)$, $B(z)F(z)=G'(z)$, i.e. $G(z)=\exp\int B(z)/A(z)dz$ and $F(z)=G(z)/A(z)$.
It remains to find/choose a suitable $C$.
Now in your case $A(z)=z^2+\lambda$ and $B(z)=z$, hence $G=\sqrt{z^2+\lambda}$ and $F=(z^2+\lambda)^{-1/2}$, i.e. we want
$$R(r)=\int_C(z^2+\lambda)^{-1/2}e^{rz}dz$$
for appropriate $C$. To get endpoints of $C$: $G(z)e^{rz}\to 0$ at $z=\pm\sqrt{-\lambda}$ and for $Re\,z\to-\infty$ (this only for $r>0$). Two natural choices for $C$ are: the segment connecting $\pm\sqrt{-\lambda}$ - this gives a solution $R(r)$ which is an entire function of $r$,  the other $C$ runs from $Re\,z=-\infty$ ($Im\, z=0$ say) to one of $\pm\sqrt{-\lambda}$, which is singular at $r=0$. ($F$ is double valued, but it only influences the sign of the integrals.) These two solutions have different behavior at $r=0$, hence they are linearly independent and form a basis of solutions. 
