In here, we have considered Bessel functions $J_\nu(x)$ for $\nu=\mbox{integer}$ case only. Note that $J_\nu$ and $J_{-\nu}$ are linearly independent if $\nu$ is a non-integer. If $\nu$ is an integer, $J_\nu$ and $J_{-\nu}$ satisfy the relation $J_{-\nu}=(-1)^\nu J_\nu$, i.e. they are no longer linearly independent. Thus we need a second solution for Bessel’s equation.

Let us define
$$N_\nu(x):=\frac{\cos\nu\pi J_\nu(x)-J_{-\nu}(x)}{\sin\nu\pi}.$$
$N_\nu(x)$ is called Neumann function or Bessel function of the second kind. For $\nu=\mbox{integer}$, $N_\nu(x)$ is an indeterminate form of type $\frac{0}{0}$. So by l’Hôpital’s rule
\begin{align*}
N_n(x)&=\lim_{\nu\to n}\frac{\frac{\partial}{\partial\nu}[\cos\nu\pi J_\nu(x)-J_{-\nu}(x)]}{\frac{\partial}{\partial\nu}\sin\nu\pi}\\
&=\frac{1}{\pi}\lim_{\nu\to n}\left[\frac{\partial J_\nu(x)}{\partial\nu}-(-1)^\nu\frac{\partial J_{-\nu}(x)}{\partial\nu}\right].
\end{align*}
Neumann function can be also written as a power series:
\begin{align*}
N_n(x)=&\frac{2}{\pi}\left[\ln\left(\frac{x}{2}\right)+\gamma-\frac{1}{2}\sum_{p=1}^n\frac{1}{p}\right]J_n(x)\\
&-\frac{1}{\pi}\sum_{r=0}^\infty(-1)^r\frac{\left(\frac{x}{2}\right)^{n+2r}}{r!(n+r)!}\sum_{p=1}^r\left[\frac{1}{p}+\frac{1}{p+n}\right]\\
&-\frac{1}{\pi}\sum_{r=0}^{n-1}\frac{(n-r-1)!}{r!}\left(\frac{x}{2}\right)^{-n+2r},
\end{align*}
where $\gamma$ is the Euler-Mascheroni number.

In Maxima, Neumann function is denoted by bessel_y(n,x). Let us plot $N_0(x)$, $N_1(x)$, $N_3(x)$ altogether using the command

plot2d([bessel_y(0,x),bessel_y(1,x),bessel_y(2,x)],[x,0.5,10]);

Neumann Functions

We now show that Neumann functions do satisfy Bessel’s differential equation. Differentiate Bessel’s equation
$$x^2\frac{d^2}{dx^2}J_{\pm\nu}(x)+x\frac{d}{dx}J_{\pm\nu}(x)+(x^2-\nu^2)J_{\pm\nu}(x)=0$$
with respect to $\nu$. (Here of course we are assuming that $\nu$ is a continuous variable.) Then we obtain the following equations:
\begin{align*}
x^2\frac{d^2}{dx^2}\frac{\partial J_\nu(x)}{\partial\nu}+x\frac{d}{dx}\frac{\partial J_\nu(x)}{\partial\nu}+(x^2-\nu^2)\frac{\partial J_\nu(x)}{\partial\nu}=2\nu J_\nu(x)\ \ \ \ \ \mbox{(1)}\\
x^2\frac{d^2}{dx^2}\frac{\partial J_{-\nu(x)}}{\partial\nu}+x\frac{d}{dx}\frac{\partial J_{-\nu(x)}}{\partial\nu}+(x^2-\nu^2)\frac{\partial J_{-\nu(x)}}{\partial\nu}=2\nu J_\nu(x)\ \ \ \ \ \mbox{(2)}
\end{align*}
Subtract $\frac{1}{\pi}(-1)^n$ times (2) from $\frac{1}{\pi}$ times (1) and then take the limit of the resulting equation as $\nu\to n$. Then we see that $N_n$ satisfies the Bessel’s differential equation
$$x^2\frac{d^2}{dx^2}N_n(x)+x\frac{d}{dx}N_n(x)+(x^2-n^2)N_n(x)=0.$$
The general solution of Bessel’s differential equation is given by
$$y_n(x)=AJ_n(x)+BN_n(x).$$

To accommodate boundary conditions for a finite interval $[0,a]$, we need to consider Bessel functions of the form $J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)$. For $x=\frac{\alpha_{\nu m}}{a}\rho$, Bessel’s equation (9) in here can be written as
\begin{equation}\label{eq:bessel10}\rho^2\frac{d^2}{d\rho^2}J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)+\frac{d}{d\rho}J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)+\left(\frac{\alpha_{\nu m}^2\rho}{a^2}-\frac{\nu^2}{\rho}\right)J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)=0.\end{equation} Changing $\alpha_{\nu m}$ to $\alpha_{\nu n}$, $J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)$ satisfies
\begin{equation}\label{eq:bessel11}\rho^2\frac{d^2}{d\rho^2}J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)+\frac{d}{d\rho}J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)+\left(\frac{\alpha_{\nu n}^2\rho}{a^2}-\frac{\nu^2}{\rho}\right)J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)=0.\end{equation}
Multiply \eqref{eq:bessel10} by $J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)$ and \eqref{eq:bessel11} by $J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)$ and subtract:
\begin{align*}
J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)\frac{d}{d\rho}&\left[\rho\frac{d}{d\rho}J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)\right]-J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)\frac{d}{d\rho}\left[\rho\frac{d}{d\rho}J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)\right]\\&=\frac{\alpha_{\nu n}^2-\alpha_{\nu m}^2}{a^2}\rho J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right).\end{align*}
Integrate this equation with respect to $\rho$ from $\rho=0$ to $\rho=a$:
\begin{equation}\begin{aligned}\int_0^\rho J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)\frac{d}{d\rho}&\left[\rho\frac{d}{d\rho}J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)\right]d\rho\\&-\int_0^\rho J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)\frac{d}{d\rho}\left[\rho\frac{d}{d\rho}J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)\right]d\rho\\&=\frac{\alpha_{\nu n}^2-\alpha_{\nu m}^2}{a^2}\int_0^\rho J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)\rho d\rho.\end{aligned}\label{eq:bessel12}\end{equation}
Using Integration by Parts, we have
\begin{align*}
\int_0^\rho &J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)\frac{d}{d\rho}\left[\rho\frac{d}{d\rho}J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)\right]d\rho\\&=\left[\rho J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)\rho\frac{d}{d\rho}J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)\right]_0^a-\int_0^a \rho\frac{d}{d\rho}J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)dJ_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right).\end{align*}
Thus \eqref{eq:bessel12} can be written as
\begin{equation}\begin{aligned}\left[\rho J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)\rho\frac{d}{d\rho}J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)\right]_0^a-\left[\rho J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)\rho\frac{d}{d\rho}J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)\right]_0^a\\=\frac{\alpha_{\nu n}^2-\alpha_{\nu m}^2}{a^2}\int_0^\rho J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)\rho d\rho.\end{aligned}\label{eq:bessel13}\end{equation}
Clearly the LHS of \eqref{eq:bessel13} vanishes at $\rho=0$. (Here we consider only $\nu=\mbox{integer}$ case.) It also vanishes at $\rho=a$ if we choose $\alpha_{\nu n}$ and $\alpha_{\nu m}$ to be $n$-th and $m$-th zeros of $J_\nu$. Therefore, for $m\ne n$
\begin{equation}\label{eq:bessel14}\int_0^\rho J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)J_\nu\left(\frac{\alpha_{\nu n}}{a}\rho\right)\rho d\rho=0.\end{equation}
\eqref{eq:bessel14} gives us orthogonality over the interval $[0,a]$.

For $m=n$, we have the normalization integral
\begin{equation}\int_0^a\left[J_\nu\left(\frac{\alpha_{\nu m}}{a}\rho\right)\right]^2\rho d\rho=\frac{a^2}{2}[J_{\nu+1}(\alpha_{\nu m})]^2.\end{equation}

Let us begin with the generating function

$$g(x,t) = e^{\frac{x}{2}\left(t-\frac{1}{t}\right)}.$$
Expanding this function in a Laurent series, we obtain
$$e^{\frac{x}{2}\left(t-\frac{1}{t}\right)} = \sum_{n=-\infty}^\infty J_n(x)t^n.$$
The coefficient of $t^n$, $J_n(x)$, is defined to be a Bessel function of the first kind of order $n$.
Now, we determine $J_n(x)$.
\begin{align*}
e^{\frac{x}{2}t}e^{-\frac{x}{2t}}&=\sum_{r=0}^\infty\left(\frac{x}{2}\right)^r\frac{t^r}{r!} \sum_{s=0}^\infty(-1)^s\left( \frac{x}{2}\right)^s \frac{t^{-s}}{s!}\\
&=\sum_{r=0}^\infty\sum_{s=0}^\infty\frac{(-1)^s}{r!s!}\left(\frac{x}{2}\right)^{r+s}t^{r-s}.
\end{align*}
Set $r=n+s$. Then for $n\ge 0$ we obtain
$$J_n(x)=\sum_{s=0}^\infty \frac{(-1)^s}{s!(n+s)!}\left(\frac{x}{2}\right)^{n+2s}.$$

Bessel Functions

Now we find $J_{-n}(x)$. In the above series for $J_n(x)$, we obtain
$$J_{-n}(x)=\sum_{s=0}^\infty\frac{(-1)^s}{s!(s-n)!} \left(\frac{x}{2}\right)^{2s-n}.$$
However, $(s-n)!\rightarrow\infty$ for $s=0,1,\cdots,(n-1)$. So the series may be considered to begin at $s=n$. Replacing $s$ by $s+n$, we obtain
$$J_{-n}(x)=\sum_{s=0}^ \infty\frac{(-1)^{s+n}}{s!(s+n)!}\left( \frac{x}{2} \right)^{n+2s}.$$ Note that $J_n(x)$ and $J_{-n}(x)$ satisfy the relation
$$J_{-n}=(-1)^nJ_n(x).$$

Let us differentiate the generating function $g(x,t)$ with respect to $t$:
\begin{align*}
\frac{\partial g(x,t)}{\partial t} &=\frac{x}{2}\left(1+ \frac{1}{t^2}\right) e^{\frac{x}{2}\left(t-\frac{1}{t}\right)}\\
&=\sum_{n=-\infty}^\infty n J_n(x) t^{n-1}.
\end{align*}
Replace $e^{\frac{x}{2}\left(t-\frac{1}{t}\right)}$ by $\sum_{n=-\infty}^\infty J_n(x) t^n$.
Then
\begin{align*}
\sum_{n=-\infty}^\infty \frac{x}{2} (1+ \frac{1}{t^2}) J_n(x) t^{n}&=\sum_{n=-\infty}^\infty \frac{x}{2} [J_n(x) t^n + J_n(x) t^{n-2}]\\
&=\sum_{n=-\infty}^\infty \frac{x}{2} [J_{n-1}(x) + J_{n+1}(x)] t^{n-1}.
\end{align*}
Thus,
$$\sum_{n=-\infty}^\infty \frac{x}{2} [J_{n-1}(x) + J_{n+1}(x)] t^{n-1}=\sum_{n=-\infty}^\infty n J_n(x) t^{n-1}$$
or we obtain the recurrence relation,
\begin{equation}J_{n-1}(x) + J_{n+1}(x) = \frac{2n}{x} J_n(x).\end{equation}
Now we differentiate $g(x,t)$ with respect to $x$:
\begin{align*}
\frac{\partial g(x,t)}{\partial x}&=\frac{1}{2}\left(1-\frac{1}{t}\right)e^{\frac{x}{2}\left(t-\frac{1}{t}\right)}\\
& =\sum_{n=-\infty}^\infty J_n'(x)t^n.
\end{align*}
This leads to the recurrence relation
\begin{equation}J_{n-1}(x) – J_{n+1}(x) = 2 J_n'(x).\end{equation}
As a special Case of this recurrence relation, we obtain,
$$J_{0}'(x)=-J_1(x).$$
Adding (1) and (2), we have
\begin{equation}J_{n-1}(x)=\frac{n}{x}J_n(x) + J_n'(x).\end{equation}
Multiplying (3) by $x^n$:
\begin{align*}
x^n J_{n-1}(x) & = n x^{n-1} J_n(x) + x^n J_n'(x)\\
& = \frac{d}{dx}[ x^n J_n(x)].
\end{align*}
Subtracting (2) from (1), we have
\begin{equation}J_{n+1}(x) = \frac{n}{x} J_n(x) – J_n'(x).\end{equation}
Multiplying (4) by $-x^{-n}$:
\begin{align*}
-x^{-n} J_{n+1}(x) & = -n x^{-n-1} J_n(x) + x^{-n} J_n'(x)\\
& = \frac{d}{dx}[x^{-n} J_n(x)].
\end{align*}

Using recurrence relations, we can show that the Bessel functions $J_n(x)$ are the solutions of the Bessel’s differential equation. The recurrence relation (3) can be written as
\begin{equation}x J_n'(x) = x J_{n-1}(x) – n J_n(x).\end{equation}
Differentiating this equation with respect to $x$, we obtain
$$J_n'(x) + x J_n^{\prime\prime}(x) = J_{n-1}(x) + x J_{n-1}'(x) – n J_n'(x)$$
or
\begin{equation}x J_n^{\prime\prime}(x) + (n+1) J_n'(x) – x J_{n-1}'(x) – J_{n-1}(x) = 0.\end{equation}
Subtracting (5) times $n$ from (6) times $x$ results the equation
\begin{equation}x^2 J_n^{\prime\prime}(x) + x J_n'(x) – n^2 J_n(x) + x(n-1) J_{n-1}(x) – x^2 J_{n-1}'(x) = 0.\end{equation}
Replace $n$ by $n-1$ in (4) and multiply the resulting equation by $x^2$ to get the equation
\begin{equation}x^2 J_n(x) = x (n-1) J_{n-1}(x) – x^2 J_{n-1}'(x).\end{equation}
With the equation (8), the equation (7) can be written as
\begin{equation}\label{eq:bessel9}x^2 J_n^{\prime\prime}(x) + x J_n'(x) + (x^2 – n^2) J_n(x) = 0.\end{equation}
This is Bessel’s equation. Hence the Bessel functions $J_n(x)$ are the solutions of Bessel’s equation.