\section{Introduction}

The need for explicit solutions of differential equations (be them ordinary or partial) arises
in the everyday practice of physicists, engineers and scientists in general. The class of ordinary
differential equations that can be solved by elementary methods, mainly clever changes of variables
and linear algebra, is quite ample and, accordingly, there exists a vast amount of software devoted
to finding both their explicit and numerical solutions. This is in contrast with the case of partial
differential equations (PDEs). Here, solution techniques are more advanced, involving Fourier analysis for
bounded domains and integral transforms in the unbounded case\footnote{In this paper we restrict our
attention to classical solutions, thus leaving aside considerations regarding weak solutions and the
use of such tools as distribution theory.} and, although there is a good deal of software packages
related to numerical computations (most  of them using the finite element method), there is a
shortage of options for finding explicit analytic solutions. In fact, to the best of our knowledge,
only Maplesoft's Maple\texttrademark\  and Wolfram's Mathematica\texttrademark\
--both closed and proprietary software-- are
capable of this task. Both are excellent general purpose computer algebra systems (CAS) and do their job
fast and efficiently, so one could wonder what is the motivation for considering the problem of
automated solution of PDEs or thinking it is not already solved. To this question we can offer the
following answers:
\begin{enumerate}
\item Although they are very good at finding solutions, still there are examples where Maple\texttrademark\
or Mathematica\texttrademark\  can not find it (see examples in subsection \ref{secex}).
In those cases, it is very difficult
to determine the origin of the failure, as both are closed source packages.
The user can have a hard time trying
to determine whether there is a problem in the formulation of the problem or in the method of solution.
\item Proprietary software usually has a cost, and this case is no exception. While both packages are worth
their price, this can be a serious hurdle for people working in underdeveloped countries, students in general
or even freelancers for whom an open source alternative would be desirable.
\item From a pedagogical point of view, the use of closed source software is like using a black box for
finding answers. The teacher can expose an algorithm in the classroom and then say it is exemplified by the
output, but really there is no clue about what is being actually computed, if there is any mathematical
trick used which is not contemplated in the algorithm, or if use is made of pre-computed tables of
particular cases, and so on.
\item Lastly, the problem is interesting in itself from a mathematical and computational perspective,
due to the blending of purely symbolic and numeric algorithms it requires, as will be seen in detail below.
\end{enumerate}

These reasons led us to write a library, called \texttt{pdefourier}, for solving the general
Sturm-Liouville problem for the heat, wave and Laplace operators with sources, as well as some particular
problems reducible to these, such as the convection-diffusion equation.
Our particular implementation uses the open source Maxima CAS \cite{Max},
which has a programming language very closely tied to LISP and should be translatable to any other CAS or
programming language in general, such as Python or Symbolic C++, without difficulty.

In Section \ref{sec2} we comment on the challenges posed by the automated computation of solutions using
Fourier techniques in general. By including examples of computations using Mathematica\texttrademark\ ,
we present a panorama of the current state of the art\footnote{We must insist in this aspect, since nothing could
be further from our minds than to do a comparison test or a benchmark against a commercial software, supported by
a whole team of scientists.}.

\begin{rem}
In what follows, all Mathematica\texttrademark\ commands have been executed using version
12.1.0.0 on a Linux platform. The Maxima commands have been executed on the same Linux machine,
with version 5.42.2 compiled against the LISP implementation SBCL 1.4.4.
It is assumed that the package \texttt{pdefourier} has already been loaded.
For installation instructions see Appendix \ref{app}.
\end{rem}

Section \ref{sec3} describes how we handle the issues mentioned in Section \ref{sec2}, detailing some of the
core features of the library, and we finish in Section \ref{sec4} by showing some examples and applications.

\section{Sturm-Liouville problems and eigenfunctions expansion}\label{sec2}

Notice that, from a mathematical point of view, the basic models given by the heat, wave and
Laplace equations are all second-order linear partial differential equations with constant
coefficients. They are supplemented by some initial and/or boundary conditions (the general situation contemplates both, and it is known as an IBVP or initial-boundary value problem), and in this
work we will assume that the problem is set on a bounded spatial domain $[0,L]\subset\mathbb{R}$,
so these conditions will be given at the boundary points $x=0,x=L$ (the only exception
being Laplace's equation, which will be considered on a plane domain).
The general problem for the heat and wave equations can be written as:
\begin{equation}\label{problem}
\begin{cases}
\mathcal{D}u =Q(x,t) \\
u(x,0)=F(x) \\
u_x(x,0)=G(x) \\
\alpha_1 u(0,t) +\beta_1 u_x(0,t)= h_1(t) \\
\alpha_2 u(L,t) +\beta_2 u_x(L,t)= h_2(t)\,,
\end{cases}
\end{equation}
where $\mathcal{D}$ denotes either the heat or the wave operator. In the case of the Laplace's
equation we consider a wider class of plane domains (rectangles, disks, wedges, annuli) and boundary conditions. In a rectangular domain \([0,a] \times [0,b] \) we consider the problem with mixed boundary conditions of the form:
\begin{equation}\label{eq-lprect}
\begin{cases}
\Delta u=0 \quad \alpha, \beta, \gamma, \delta \in \{ 0,1 \} \\
(1-\alpha)u(x,0)+\alpha u_y(x,0)=f_0(x) \\
(1-\beta)u(x,b)+\beta u_y(x,b)=f_b(x) \\
(1-\gamma) u(0,y) + \gamma u_x (0,y)=g_0(y) \\
(1-\delta) u(a,y) + \delta u_x (a,y)= g_a(y)
\end{cases}
\end{equation}

On the remaining domains, the Laplace equation is solved in polar coordinates with Dirichlet and Neumann boundary conditions; the only exception being the case of an annular domain, where only Dirichlet boundary conditions are considered. Here is a
summary of the cases:

\begin{equation}\label{eq-lpdwa}
\begin{cases}
u_{rr}+\frac{1}{r}u_r+\frac{1}{r^2}u_{\theta\theta}=0 \\
(r,\theta) \in (0,R_1]\times [0,\alpha[, \, 0<\alpha\leq 2\pi, \,0<R_1 \text{ (Disk/Wedge)} \\
(r,\theta) \in [R_1,R_2]\times [0,2\pi[, \, 0<R_1<R_2 \text{ (Annulus)} \\
\text{Dirichlet boundary conditions:} \\
u(r,0)=0=u(r,\alpha) \text{ if } \alpha<2\pi \text{ (Wedge)} \\
u(R_1,\theta)=f(\theta)  \text{ (Annulus/Disk/Wedge)}\\
u(R_2,\theta)=g(\theta) \text{ (Annulus)} \\
\text{Neumann boundary conditions:} \\
u(r,0)=0=u(r,\alpha) \text{ if } \alpha<2\pi \text{ (Wedge)}\\
u_r(R_1,\theta)=f(\theta) \text{ (Disk/Wedge)}\\
\end{cases}
\end{equation}

A common procedure to solve these equations is to apply separation of variables: Writing
$$
u(x,t)=X(x)T(t)\,,
$$
and denoting generically $X$ or $T$ by $\phi$, leads to a regular Sturm-Liouville problem
whose general form is
\begin{equation}\label{SL}
\mathcal{L}\phi \doteq (p\cdot\phi')'+(q +\lambda r)\phi =0\,,
\end{equation}
where $p,r>0$ are positive real functions, and $p',q,r$ are continuous on $[0,L]$. The abstract theory of Sturm-Liouville operators is very well developed (see \cite{AG08,Hel67}), and its main result can be summarized as follows:

\begin{thm}
Let $\mathcal{L}$ be a Sturm-Liouville operator as in \eqref{SL}, subject to the boundary conditions
derived from those in \eqref{problem} in the homogeneous case. Assume that, for $i\in\{1,2\}$,
\begin{equation}\label{conds}
|\alpha_i|^2+|\beta_i|^2\neq 0\,.
\end{equation}
Then, the following hold:
\begin{enumerate}
\item All the eigenvalues $\lambda$ are real and form a denumerable set.
\item The set of all eigenfunctions, $\mathcal{E}$, is an orthogonal set.
\item The set $\mathcal{E}$ is complete in the space of square-integrable functions
$L^2([0,L])$.
\item The sequence $(\lambda_n)$ diverges to $+\infty$ and is bounded from below.
\item If $\alpha_i,\beta_i,p,-q$ are all non-negative, then the eigenvalues $\lambda_n$ are
all non-negative.
\end{enumerate}
\end{thm}

Notice that the preceding theorem applies to \emph{homogeneous} initial and boundary conditions and \(Q=0 \).
To deal with the non-homogeneous case , we use the method of eigenfunction expansion or finite cosine or sine transform (see \cite{Hab13,Wal98}).
Once we know how to deal with the general Sturm-Liouville problem for the heat, wave and Laplace
operators, we can solve some other problems that can be reduced to one of these. The cases of the
convection-difussion equation and heat loss through lateral boundaries are then particular cases.
Thus, consider the general expression of a second-order linear parabolic equation with constant
coefficients with sources:
\begin{equation}\label{auxiliar}
u_t=\kappa u_{xx} + v u_x + cu +Q(x,t)
\end{equation}
Let us apply the change of variable
\[
u=we^{\rho(x,t)}\,,
\]
where
\[
\rho(x,t)=(c-\frac{v^2}{4\kappa})t-\frac{v}{2\kappa}x
\]
Substituting into \eqref{auxiliar} we get
\begin{align*}
Q&=u_t-\kappa u_{xx} -vu_x-cu \\
&=e^\rho(w_t-\kappa w_{xx})
\end{align*}
hence, the original equation in terms of the new dependent variable reduces to the heat equation with sources
\[
w_t=\kappa w_{xx}+e^{-\rho}Q
\]

Of course, the initial and boundary conditions for $u(x,t)$ must be transformed accordingly into conditions
for $w(x,t)$.

The Sturm-Liouville theory is the basic tool used in our Maxima package \texttt{pdefourier}. The algorithm
can be summarized as follows:  we first
solve the eigenvalue problem associated to \eqref{SL}, normalize the corresponding eigenfunctions and then proceed to develop the solution to \eqref{problem} in term of eigenfunctions, taking non-homogeneous conditions into account with the use of the methods commented above. These steps involve some delicate issues that we discuss next.

\textbf{Singular values of the Fourier coefficients}. The computation of Fourier coefficients in a CAS can be done using integration formulas. Particular care must be taken to separate singular values that lead to divergences in the general expressions of the Fourier coefficients. There are two types of special input functions whose Fourier coefficients must be calculated separately, namely, the functions belonging to one of the cosine, sine, trigonometric or complex orthogonal systems or a product of a polynomial with them. In these cases, the Fourier coefficients have a general expression
$$
\frac{C}{\rho (n)}
$$
with $\rho$ a polynomial, which is only valid for the non-singular values of $n$. The singular
values must be treated separately. This process is, of course, transparent to the user, but lies at the bottom of the difficulties experienced by such mature software as Mathematica\texttrademark\ . Due to the
widespread of this CAS, we will use it to get a benchmark for our library performance.


To get a sense of the kind of problems we are talking about, recall the definition of the Fourier
coefficients of a $2\pi-$periodic function $f(x)$ (we consider the case $L=\pi$ for simplicity):

\[
a_0=\frac{1}{2\pi}\int^\pi_{-\pi}f(x)\,\mathrm{d}x\,,\
a_n=\frac{1}{\pi}\int^\pi_{-\pi}f(x)\cos (nx)\,\mathrm{d}x\,,\
b_n=\frac{1}{\pi}\int^\pi_{-\pi}f(x)\sin (nx)\,\mathrm{d}x\,.
\]
where $n$ is a natural number.
In practice, the explicit computations of these coefficients depends heavily on the orthogonality
relations ($m,n\geq 1$)
\[
\frac{1}{\pi}\int^\pi_{-\pi}\cos (mx)\cos (nx)\,\mathrm{d}x =\delta_{mn}\,,\
\frac{1}{\pi}\int^\pi_{-\pi}\sin (mx)\sin (nx)\,\mathrm{d}x =\delta_{mn}\,.
\]

However, most computer algebra systems fail to recognize the possibility that $m=n$. Thus,
Mathematica\texttrademark\ gives
\begin{verbatim}
In[1]:= Integrate[Cos[m x] Cos[n x], {x, -Pi, Pi}]
Out[1]= (2 m Cos[n \[Pi]] Sin[m \[Pi]] - 2 n Cos[m \[Pi]] Sin[n \[Pi]])/(m^2 - n^2)
\end{verbatim}
and the same is true of Maxima:
\begin{verbatim}


\end{verbatim}

This fact has many consequences due to the presence of singular values of $n$, as pointed above.
Take as an example the computation of Fourier coefficients for
\begin{equation}\label{fexample}
f(x)=3x^2\cos(7x)\,.
\end{equation}
Using the trigonometric identity
\[
\cos\theta\cos\phi =\frac{1}{2}\cos(\theta-\phi)+\frac{1}{2}\cos(\theta +\phi)
\]
and integrating by parts, it is straightforward to obtain
\begin{align*}
a_0=& -\frac{6}{49} \\
a_n=& \begin{cases}
-\dfrac{12(n^2+49)(-1)^n}{(n-7)^2(n+7)^2}\mbox{ if }n\neq 7 \\[7pt]
\dfrac{98\pi^2+3}{98}\mbox{ if }n=7
\end{cases}\\
b_n=& 0
\end{align*}

Mathematica\texttrademark\  fails to compute correctly the Fourier coefficients directly from the integral formulas because it does not recognize the singularity corresponding to $n=7$:
\begin{verbatim}
In[2]:= Assuming[Element[n,Integers],1/Pi Integrate[Cos[n x] 3 x^2 Cos[7 x],{x,-Pi,Pi}]]
Out[2]= -((12 (-1)^n (49 + n^2))/(-49 + n^2)^2)
\end{verbatim}
and so does Maxima:
\begin{verbatim}

        factor(integrate(cos(n*x)*3*x^2*cos(7*x),x,

\end{verbatim}

\textbf{Equivalence of trigonometric expressions}.
Another source of concern when dealing with trigonometric expressions is the existence of multiple, equivalent ways of writing
them.
This is a well-known issue in different CASs.
The set of heuristics rules Mathematica\texttrademark\, appears to use inside its Fourier coefficients built-in functions try to support common special cases, but are highly sensitive to the way the input is written and do not take advantage of the linearity of the integral nor distribute over subintervals of piecewise-defined functions. For instance, using equation \eqref{fexample} in the case $\theta=\varphi$ we get:
\begin{equation}\label{trigexample}
\cos (\theta)^2 = \frac{1}{2}+\frac{1}{2}\cos{\left( 2 \theta \right) }
\end{equation}

For illustration purposes, here are the results given by Mathematica\textsuperscript{TM} for the Fourier cosine coefficients of each side of the equation.
\begin{verbatim}
In[3]:=FourierCosCoefficient[Cos[x]^2, x, n]
Out[3]:=0
In[4]:=FourierCosCoefficient[(1+Cos[2x])/2, x, n]
Out[4]:=1/2 (DiscreteDelta[-2 + n] + 2 DiscreteDelta[n])
\end{verbatim}

Moreover, its routines do not take advantage of the linearity, as a slight change in the input function prevents them from
working:
\begin{verbatim}
In[3]:=FourierCosCoefficient[(2+Cos[2 x])/2,x,n]
Out[3]:=0
\end{verbatim}

The strategy followed in our implementation of the package is to internally transform any trigonometric
function into its canonical form, so it becomes easy to decide whether or not the input contains an
expression whose Fourier coefficients have singular values using the pattern matching capabilities of the
Maxima CAS, as described in the following section.

\textbf{Piecewise-defined functions}. Let us point that, in applications, many of the functions we deal with are
piecewise defined, and this poses its own challenges. Maxima currently does not have built-in support for piecewise-defined expressions which are of particular importance in engineering and physics applications. We developed another package called \texttt{piecewise} automatically loaded by \texttt{pdefourier} to deal with this kind of functions.  In particular, we are able to detect the parity of a piecewise-defined function defined on an interval \([-L,L] \) by the means of comparison of subintervals allowing us to simplify the computation of Fourier coefficients of odd or even functions. Also, we take advantage of the linearity and interval addition properties of the integral; because of this, our pattern matching rules work even with piecewise-defined functions and equivalent trigonometric expressions as we will see in Section \ref{sec4}. \\



\section{Computation of the Fourier coefficients and series}\label{sec3}

The aim of this section is to explain the strategies followed in the implementation of our package to tackle the different challenges we have showed in the previous section.

\subsection{Piecewise-defined functions}

Inherited from Lisp, Maxima's main data structures are lists, thus the natural way to proceed was to translate a piecewise-defined function to a list and viceversa. Having a way to convert a piecewise expression to a list, it is straightforward to write some other functions to make operations between them or to compute the derivative and the integral of a piecewise expression written as a list. Although the package \texttt{piecewise} is able to work in more general settings, we restrict ourselves here to show only how the package works in a bounded interval. The remaining details are available in the package
documentation. For illustration purposes, consider an if-else expression of the form: \\

\texttt{if x>=a$_0$ and x<=a$_1$ then expr$_1$ elseif ... elseif x>a$_{n-1}$ and x<=a$_{n}$ then expr$_{n}$ } \\

We impose the restrictions $a_i\leq a_{i+1}, a_i \in \mathbb{R},$ to avoid potential issues with the detection of valid intervals in all the equivalent ways of writing them, and to do it efficiently.
The output returned by \texttt{pw2list} has the form:
\[
[[[a_0,a_1],\text{expr}_1],[[a_1,a_2],\text{expr}_2,]\ldots]
\]

The approach of internally working with lists instead of if-else expressions is particularly handy since the Maxima
language supports the functional programming paradigm so once we have proper detection of special cases with pattern-matching,
we can easily map the routines to each subinterval.

\subsection{Computation of Fourier coefficients}


As mentioned in the preceding Section, the product of a polynomial by a trigonometric function gives rise to
singular values when determining Fourier coefficients, having as particular cases trigonometric functions whose
wave number is a multiple of $\pi/L$. Because of the linearity of the integral, using trigonometric canonical forms
it is enough to detect the following patterns:
\[x^r\cos (\frac{m \pi x}{L}), x^r\sin (\frac{m \pi x}{L}), \cos (\frac{m \pi x}{L}),\sin (\frac{m \pi x}{L}), \quad r,m \in \mathbb{N}\]



To do so, we used the Maxima built-in commands \texttt{defmatch} and \texttt{matchdeclare}. Then, two similar strategies were followed depending on the input expression being piecewise-defined or not. Algorithm \ref{algorithm:fcoeffs} describes the steps followed in our implementation to calculate the Fourier coefficients.
Of course, when the input function is piecewise-defined the procedure described in the algorithm must be applied to each subinterval of the domain.

\begin{algorithm}[htp]\label{algorithm:fcoeffs}
\SetAlgoLined
\DontPrintSemicolon
\KwInput{ expression \emph{expr}; variable \emph{var}, semi-length of interval \emph{L} }
\KwOutput{ A list of Fourier coefficients of the form \([[a_0,a_n,b_n], l.s.v ]\) \\ l.s.v=list of singular values}
Apply simplification functions to convert trigonometric expressions appearing in \emph{expr} into their canonical form.\;
Expand \emph{expr} fully and convert the expanded expression \( term_1+\ldots+term_s\) into a list \([term_1,\ldots,term_s ] \). \;
Apply a heuristic routine with pattern matching capabilities to each element of the list to compute \( a_n(term_k), \, b_n(term_k) \) and apply an auxiliary function which searches the set of indices having singular values, if any, of \( a_n, b_n \) and store it in a list \(A\).\;
Sum over \( k\) both \( a_n(term_k), \, b_n(term_k) \) to obtain the final answer for \(a_n, b_n \) and compute separately the Fourier coefficients of the indices appearing in \(A\) storing the result in the list of singular values. The list of singular values, if any, is written as follows \( [[j,a_j,b_j],\ldots] \), otherwise, an empty list [\,] is returned. \;
\caption{Computation of Fourier coefficients}
\end{algorithm}


The same ideas also apply for the case of the complex, sine or cosine coefficients; the only difference is in the way the output is returned to the user (see Table \ref{table:formats}).

\begin{table}[pos=h!]
\centering
\begin{tabular}{llll}
Type & Command & Answer format & List of s.v format \\
\hline\\\\[-1.9\medskipamount]
Trigonometric & \texttt{fouriercoeff(expr, var, L)} & $[[a_0, a_n, b_n], \mbox{list of s.v}]$ & $[\,] \mbox{ / } [[j,a_j,b_j],\ldots]$ \\
Complex &  \texttt{cfouriercoeff(expr, var, L)} & $[[c_0,c_n], \mbox{list of s.v}]$ & $[\,] \mbox{ / } [[j,c_j],\ldots]$ \\
Cosine & \texttt{fouriercoscoeff(expr, var, L)} & $[[a_0,a_n], \mbox{list of s.v}]$ & $[\,] \mbox{ / } [[j,a_j],\ldots]$ \\
Sine & \texttt{fouriersincoeff(expr, var, L)} & $[[b_n], \mbox{list of s.v}]$ & $[\,] \mbox{ / } [[j,b_j],\ldots]$ \\[7pt]
\hline
\end{tabular}
\caption{Output formats for the Fourier coefficients (s.v= singular values).}
\label{table:formats}
\end{table}


Also notice that, by using expressions instead of declared functions, we are capable to treat arbitrary,
non-evaluated functions such as a symbolic $f(x)$. The capabilities of Maxima, however,
allows us to evaluate these when required.

\subsection{Fourier series}

Following our policy of efficiency, Fourier series are obtained using an expansion routine of the Fourier coefficients. We can suspect that Mathematica\texttrademark\, does not perform an expansion of the coefficients to obtain the Fourier series, but computes them one by one inside the \texttt{FourierSeries} function, for two reasons: First, because Mathematica\texttrademark\, only works with truncated series and so a general expression for the coefficients is not needed and second, because this approach avoids the problems associated to simplification of symbolic integrals using assumptions. Nevertheless this strategy is inefficient because it does not take advantage of previously computed Fourier coefficients. As we have seen before, Mathematica\texttrademark\, can not obtain correctly the cosine coefficients of \( \cos^2 x \), but if we compute the first five terms of the cosine series of this function we get the right answer:
\begin{verbatim}
In[10]:=FourierCosSeries[Cos[x]^2,x,5]
Out[10]:=1/2+1/2Cos[2x]
\end{verbatim}

In our case the upper limit of summation can be a positive integer or infinite. In the first case, a truncated series is returned; in the second, a symbolic series is displayed. In Table \ref{table:expansion} we summarize the syntax for the expansion routines along with the different Fourier series commands, see also Table \ref{table:formats} for comparison.
\begin{table}[pos=htp]
\centering
\begin{tabular}{l|l}
Expansion routine & Series \\
\hline\\\\[-1.9\medskipamount]
\texttt{fouriercoeff\_expand(list of coeff,var,L,N)} & \texttt{fourier\_series(expr,var,L,N)} \\
\texttt{cfouriercoeff\_expand(list of coeff,var,L,N)} & \texttt{cfourier\_series(expr,var,L,N)} \\
\texttt{fouriersincoeff\_expand(list of coeff,var,L,N)} & \texttt{fouriersin\_series(expr,var,L,N)} \\
\texttt{fouriercoscoeff\_expand(list of coeff,var,L,N)} & \texttt{fouriercos\_series(expr,var,L,N)} \\[7pt]
\hline
\end{tabular}
\caption{Expansion routines and Fourier series syntax in \texttt{pdefourier}.}
\label{table:expansion}
\end{table}

\subsection{Numerical aspects}
The determination of Fourier coefficients in the case of the $2D$ wave equation requires computing the zeros of Bessel's
functions. Unfortunately, Maxima 's numerical capabilities are rather limited, in particular, it does not have any built-in
function for this task (although it has implemented the Bessel functions of first and second kind). Thus, it has been
necessary to write a dedicated function for computing these zeros, and in the process of doing so several surprises
regarding the implementation of Bessel functions in commercial CASs arose.

Our solution for determining the zeros of $J_\nu(x)$ and $J'_\nu(x)$ is based on the papers \cite{GZ73} and \cite{IK91},
which allows for a simple and efficient implementation of their algorithms, the main step being the computation of
eigenvalues of certain matrices (we do this by importing the function \texttt{dgeev} from Netlib's \texttt{lapack}
package). As we are interested exclusively in real zeros, we restrict the order in $J_\nu(x)$ to $\nu>-1$ (Lommel's theorem),
but otherwise we impose no restrictions on $\nu$. Thus, for example, the third zero of $J_{1/2}(x)$ is returned by
the function \texttt{BesselJzeros} (mimicking the names in Maple\texttrademark\, and Mathematica\texttrademark\,):
\begin{verbatim}



\end{verbatim}
and the third zero of $J_{146225}(x)$ is
\begin{verbatim}


\end{verbatim}

While Maple\texttrademark\, gives the result (even with better precision) effortlessly, to our surprise
Mathematica\texttrademark\, could not calculate it, neither did the web-based engine Wolfram Alpha.
Probably, this is due to the fact that they are using some variant
of MacMahon formula, which fails for higher orders (see the comments in \cite{Wat18}).
Other alternatives, such as Halley's method, also put a bound on the order, as it occurs in Matlab \cite{Nic20}.

As mentioned, the algorithms in \cite{GZ73,IK91} can be used to generate the zeros of $J'_\nu(x)$. This is done
by the function
\texttt{BesselJdiffZeros}. The following table gives the first 5 zeros of the derivatives of $J'_\nu(x)$ for
$0\leq\nu\leq 4$, and is to be compared with the one in \cite{Wei20}:
\begin{verbatim}


		[3.831705970207512,	1.84118378134066,	3.05423692822714,	4.201188941210528,	5.317553126083997],
		[7.015586669815619,	5.331442773525031,	6.706133194158456,	8.015236598375953,	9.28239628524161],
		[10.17346813506272,	8.536316366346284,	9.96946782308759,	11.345924310743,	12.68190844263889],
		[13.32369193631421,	11.70600490259207,	13.17037085601612,	14.58584828616704,	15.96410703773154],
		[16.47063005087759,	14.86358863390901,	16.34752231832178,	17.78874786606648,	19.19602880004888]
	)
\end{verbatim}

Up to our knowledge, Maple\texttrademark\, does not have a similar function for computing the zeros of derivatives of
Bessel functions.


\subsection{Solution of PDEs}
In developing our library, we focused on solving the three main second-order linear partial differential equations with constant coefficients, namely, the heat, wave and Laplace equations.
The mathematical details can be found in any standard textbook on PDEs or Fourier Analysis (for instance,
see \cite{Hab13,Wal98}). Here, we just briefly discuss some technical aspects of the implementation, a list of examples will be given in the next Section.

It has been mentioned that Fourier series are obtained by performing an expansion of the Fourier coefficients. This approach is multipurpose. It is also useful to code the solution of the PDEs in an easier way. Although the package can solve the heat equation with a heat source $Q(x,t)$, we want to illustrate how the expansion sub-routines facilitate the process of finding solutions with a simpler instance of the general equation. Consider the following IBVP:
\[
\begin{cases}
u_t = \kappa u_{xx} \qquad (x,t) \in [0,L] \times \mathbb{R^+} \\
u(x,0)=f(x) \\
u(0,t)=0 \\
u(L,t)=0
\end{cases}
\]

the solution is then given by
\begin{align*}
u(x,t) &= \sum_{n=1}^ \infty B_n \exp [-n^2 ( \frac{\pi \kappa}{L})^2 t ] \sin \frac{n \pi x}{L} \\
B_n&=\frac{2}{L} \int_0^L f(x) \sin \frac{n \pi x }{L}
\end{align*}

Notice that for our purposes, it will be sufficient if we create a list of the form $[[C_n(t)],[ \, ] / [j,C_j(t)] ]$ where $C_n(t)=B_n \exp [-n^2 ( \frac{\pi \kappa}{L})^2 t ]$ and then use the expansion routine corresponding to a Fourier sine series on the interval $[0,L]$ and the space variable $x$. For instance, if we consider $f(x)=x^2(1-x),  \, L=1, \, \kappa=1$, a solution to the IBVP can be obtained in a few lines of code:

\begin{verbatim}





\end{verbatim}
\[
-\frac{4 \sum_{n=1}^{\infty }{\left. \frac{\left( 2 {{\left( -1\right) }^{n}}+1\right) \, {{\% e}^{-{{\ensuremath{\pi} }^{2}}\, {{n}^{2}} t}} \sin{\left( \ensuremath{\pi}  n x\right) }}{{{n}^{3}}}\right.}}{{{\ensuremath{\pi} }^{3}}}
\]

Here the list of singular values was empty (obviously), and no further work was required. The library
command \texttt{mixed\_heat} automates all the process in the most general setting.


Similar methods to obtain solutions for the three equations with different types of boundary conditions and domains have been implemented, dealing with general expressions (piecewise-defined or not) and taking care of the possible singular values in the coefficients.


\section{Some examples}\label{sec4}

\subsection{Fourier coefficients and series}

Let $f(x)$ be as in \eqref{fexample}. Then we have the following:

\begin{verbatim}


\end{verbatim}
\[
[[-\frac{6}{49},\frac{12 \left( {{n}^{2}}+49\right) \, {{\left( -1\right) }^{n+1}}}{{{n}^{4}}-98 {{n}^{2}}+2401},0],[[7,\frac{98 {{\ensuremath{\pi} }^{2}}+3}{98},\frac{2}{7}]]]
\]

Indeed, notice that $(n-7)^2(n+7)^2=n^4-98n^2+2401$.

The case where Mathematica\texttrademark\,  fails to compute correctly the cosine coefficients of equivalent trigonometric expressions \eqref{trigexample}, is readily solved by our package:

\begin{verbatim}


\end{verbatim}
\[
[[\frac{1}{2},0],[[2,\frac{1}{2}]]]
\]

\begin{verbatim}


\end{verbatim}
\[
[[\frac{1}{2},0],[[2,\frac{1}{2}]]]
\]

Next, we offer an example of a piecewise-defined function having singular values. We will use the command \texttt{cfouriercoeff} to get the complex Fourier coefficients, so that we can compare our answer to the one returned by Mathematica\texttrademark.

\begin{verbatim}



\end{verbatim}
\[[[\frac{1}{3 \ensuremath{\pi} },-\frac{3 \left( {{\left( -1\right) }^{n}}+1\right) }{2 \ensuremath{\pi}  \left( {{n}^{2}}-9\right) }],[[3,-\frac{\% i}{4}]]]\]

However, Mathematica \textsuperscript{TM} is not able to detect the singular value of the coefficient when $n=3$:
\begin{verbatim}
In[8]:=f[x]:= Piecewise[{{0,-Pi<=x<0}, {Sin[3x],0<=x<Pi}}]
In[9]:=FourierCoefficient[f[x], x, n]
\end{verbatim}
\[
-\frac{3 \left((-1)^n+1\right)}{2 \pi  \left(n^2-9\right)}
\]

Now, we consider an example of how Fourier series are displayed symbolically. If $g(x)=x^4$ on $[-\pi,\pi]$ we get the answer in the form used in textbooks:

\begin{verbatim}


\end{verbatim}
\[
8 \left( \sum_{n=1}^{\infty }{\left. \frac{\left( {{\ensuremath{\pi} }^{2}}\, {{n}^{2}}-6\right) \, {{\left( -1\right) }^{n}} \cos{\left( n x\right) }}{{{n}^{4}}}\right.}\right) +\frac{{{\ensuremath{\pi} }^{4}}}{5}
\]

Lastly, we show an example of how to use the expansion routines to obtain the Fourier series, truncated or not, and how they handle singular values in the Fourier coefficients when displaying an infinite series.

\begin{verbatim}

(fcoeff)
\end{verbatim}
\[
[[\frac{2 n\, {{\left( -1\right) }^{n}}}{{{n}^{2}}-9}],[[3,-\frac{1}{6}]]]
\]
\begin{verbatim}


\end{verbatim}
\[
-\frac{5 \sin{\left( 5 x\right) }}{8}+\frac{8 \sin{\left( 4 x\right) }}{7}-\frac{\sin{\left( 3 x\right) }}{6}-\frac{4 \sin{\left( 2 x\right) }}{5}+\frac{\sin{(x)}}{4}
\]
\begin{verbatim}


\end{verbatim}
\[
2 \left( \sum_{n=1}^{\infty }{\left. \frac{n\, {{\left( -1\right) }^{n}} \sin{\left( n x\right) }}{{{n}^{2}}-9}\right.}\right) -\frac{\sin{\left( 3 x\right) }}{6}
\]

It is important to notice that displaying infinite series correctly has been a source of troubles in different CAS.
This is mainly due to the fact that they are not evaluated, only displayed symbolically, so
a simplification evident for a human might not be performed by a CAS.

For illustration purposes, suppose that we want to compute the Fourier series of $h(x)=\sin 15 x$ on the interval $[-\pi,\pi]$.
It is obvious that its Fourier series is exactly equal to $h(x)$, because its Fourier coefficients are $a_0=0, \, a_n=0, \, b_n=\delta_n^{15}$ and we have:
\[
h(x) \sim \sum_{n=1}^\infty \delta_n^{15} \sin n x = \sin 15 x
\]

However, infinite sums of expressions containing Kronecker delta functions are really hard to simplify to a
single term, because it must be verified that the only index that does not vanish is indeed contained in the
set of indices over which we are considering the sum. Integration routines in different CAS sometimes return
the result for \(b_n \) in terms of some sort of Kronecker delta (Mathematica \textsuperscript{TM} returns it in
terms of \texttt{DiscreteDelta}), which complicates the task of obtaining the Fourier series by an expansion
of the coefficients. Following our approach, we avoid the issue of the evaluation of Kronecker delta functions inside an infinite sum:

\begin{verbatim}


\end{verbatim}

Notice that, as we just saw in the examples above, when the list of singular values is not empty, the
library prints a message warning about the set of indices excluded in the infinite sum.

\subsection{Partial differential equations}\label{secex}

We will show some examples involving each type of equation (parabolic, hyperbolic and elliptic).
Although \texttt{pdefourier} is able to solve these equations with many more boundary conditions, the ones
given here are meant to give a panorama about how to use the package to solve PDEs.

\subsubsection{The linear second-order parabolic equation}
Consider the general problem:
\[
\begin{cases}
u_t=\kappa u_{xx} + v u_x + cu +Q(x,t) \mbox{ on } (x,t) \in [0,L] \times \mathbb{R}^+ \\
u(x,0)=F(x) \\
\alpha_1 u(0,t)+\beta_1 u_x(0,t) =h_1(t) \\
\alpha_2 u(L,t)+\beta_2 u_x(L,t) =h_2(t)
\end{cases}
\]

The syntax for solving it is
\[
\mbox{\texttt{mixed\_parabolic(Q(x,t),F(x),a1,b1,a2,b2,h1(t),h2(t),x,t,L,k,v,c,ord)}},
\]
where \texttt{ord} can be \texttt{inf}, for the complete series solution, or a natural number, for
the truncated series.

As an example, we will solve a problem with homogeneous boundary conditions
\[
\begin{cases}
u_t=\kappa u_{xx} - 9 u_x \mbox{ on } (x,t) \in [0,1] \times \mathbb{R}^+ \\
u(x,0)=e^{45x/10}(5\sin(\pi x)+9\sin(2\pi x)+2\sin(3\pi x)) \\
u(0,t) = 0) \\
u(1,t) =0
\end{cases}
\]
\begin{rem}
In comparing results with Maple\texttrademark\, 2020 and Mathematica\texttrademark\, we
imposed a time limit of $600$ seconds (real time) to get an answer from both. If no answer was returned in such time,
we considered it a failure. Maple\texttrademark\, reached that time limit without giving an answer;
however, Mathematica\texttrademark\, solved it giving the same answer as below.
\end{rem}

With \texttt{pdefourier} we would do
\begin{verbatim}


                        1,0,1,0,0,0,x,t,1,1,-9,0,inf);

\end{verbatim}
\[{e^{\frac{9 x}{2}-\frac{81 t}{4}}}\, \left( 2 {e^{-9 {{\ensuremath{\pi} }^{2}} t}} \sin{\left( 3 \ensuremath{\pi}  x\right) }+9 {e^{-4 {{\ensuremath{\pi} }^{2}} t}} \sin{\left( 2 \ensuremath{\pi}  x\right) }+5 {e^{-{{\ensuremath{\pi} }^{2}} t}} \sin{\left( \ensuremath{\pi}  x\right) }\right) \]



\subsubsection{Laplace equation}
Here we consider the Laplace equation on a wedge with Neumann conditions:
\[
\begin{cases}
u_{rr}+\frac{1}{r} u_r + \frac{1}{r^2}u_{\theta \theta}=0 \mbox{ on } (r,\theta) \in [0,R] \times [0,\alpha] , \, 0<\alpha < 2\pi \\
u(r,0)=0 \\
u(r,\alpha)=0 \\
u_r(R,\theta)=f(\theta)
\end{cases}
\]

Mathematica\texttrademark\, is unable of solving it (tested on versions 12.0 and 12.1).
The syntax for Neumann problems in \texttt{pdefourier} is
\[
\mbox{\texttt{neumann\_laplace\_wedge(R,alpha,f(theta),theta,ord)}}.
\]

In this example we will consider\footnote{The use of the percent sign is useful in some
Maxima's graphical interfaces to get the symbol displayed as a Greek character.} $\alpha= \pi/2,R=1$, and an arbitrary
function $f(\theta)$:

\begin{verbatim}


\end{verbatim}
\[
\frac{2}{\pi} \sum_{n=1}^{\infty }{\left. \frac{\sin{\left( 2 \theta  n\right) } \int_{0}^{\frac{\ensuremath{\pi} }{2}}{\left. \operatorname{f}\left( \theta \right)  \sin{\left( 2 \theta  n\right) }d\theta \right.}\, {{r}^{2 n}}}{n}\right.}
\]

\subsubsection{Wave equation}
For this example, we will consider the wave equation with homogeneous boundary conditions and a driving term,
a problem whose general expression is:
\[
\begin{cases}
u_{tt}(x,t) = c^2 u_{xx}(x,t)+F(x,t) \quad (x,t) \in [0,L] \times \mathbb{R}^+\\
u(x,0)=f(x) \\
u_t(x,0)=g(x) \\
u(0,t)=0 \quad u(L,t)=0
\end{cases}
\]

The syntax for solving it using \texttt{pdefourier} is
\[
\mbox{\texttt{dirichlet\_wave(F(x,t),f(x),g(x),h\_1(t),h\_2(t),x,t,L,c,ord)}}.
\]

If $F(x,t)= \cos(\omega t) r(x)$ and $g(x)=0$, then:

\begin{verbatim}



\end{verbatim}
\begin{multline}
\notag
- \frac{\sum_{n=1}^\infty \sin ( \frac{\pi n x}{L }) \left( (2L^2 \cos ( \frac{\pi c n t }{L}) - 2L^2 \cos (\omega t) ) \int_0^L r(x) \sin (\frac{\pi n x}{L}) dx \right)}{\pi^2Lc^2 n^2 -L^3 \omega^2}  \\
-\frac{ \sum_{n=1}^\infty ( 2L^2 \omega^2-2\pi^2c^2 n^2)\sin ( \frac{\pi n x}{L }) \cos ( \frac{\pi c n t }{L}) \int_0^L f(x) \sin (\frac{\pi n x}{L}) dx}{\pi^2Lc^2 n^2 -L^3 \omega^2}
\end{multline}

Since no assumptions were made about $\omega$, the solution corresponds to the case without resonance ($\omega \neq c n \pi/L$). This is Exercise 8.5.2(b) in \cite{Hab13}. For a different example, consider a vibrating clamped circular membrane,
which is modeled in polar coordinates as
\[
\begin{cases}
u_{tt}=c^2\left(u_{rr}+\frac{1}{r} u_r +\frac{1}{r^2}u_{\theta \theta}\right)\mbox{ on }(r,\theta,t)\in [0,R]\times [0,2\pi]
\times [0,\infty [  \\
u(r,\theta,0)=f(r,\theta) \\
u_t(r,\theta,0)=g(r,\theta) \\
u(a,\theta,t)=0 \\
u(r,0,t)=u(r,2\pi,t) \\
u_\theta (r,0,t)=u_\theta (r,2\pi,t) \
\end{cases}
\]

The general case is handled by the function
\[
\mbox{\texttt{wave2d\_disk(c,R,f,g,k,l)},}
\]
where $k$ is the maximum order of the Bessel functions desired in the solution and $l$ is the number of positive zeros
taken into account.
In particular, for a circular membrane of unit radius, $c=1$, initial shape described by $f(r,\theta)=1-r^4$,
and no initial velocities, we would issue the following commands:
\begin{verbatim}



\end{verbatim}

\begin{rem}
There is a difference between this function and the remaining ones in the library. While we used expressions
as arguments before (as \texttt{sin(5*\%pi*x)}), here we need to enter declared functions (like \texttt{f,g}).
Ultimately, this has to do with the fact that we are using numerical computations, but this is not an essential
aspect, and as such will surely be modified in future versions to allow for expressions in the arguments.
\end{rem}

We have suppressed the output of the last command, since it contains a lot of terms with the aspect
\[
J_j(3,6{.}380161895923982\,r)\cos(6{.}380161895923982\,t) 2{.}467777882600071\cdot 10^{-15}\sin(3\theta)\,.
\]
Obviously, these are numerical artifacts coming from rounding errors (as revealed by exponents like $10^{-15}$).
In these cases, we can employ the function \texttt{chop} (provided in the package) analogous to
	the \texttt{Chop[]} encountered in Mathematica\texttrademark\,, to avoid rounding errors in the plot function.
	By default, \texttt{chop} discards terms of absolute value
	less than $10^{-12}$. An additional argument can be supplied, as in \texttt{chop(expr,exponent)}, to
	chop those terms whose absolute value is less that $10^{-exponent}$:
\begin{verbatim}

(expr)	1.366663216985716*bessel_j(0.0,2.404825557695772*r)*cos(2.404825557695772*t)
-0.4858370155994775*bessel_j(0.0,5.520078110286311*r)*cos(5.520078110286311*t)
\end{verbatim}

The output has a format suitable for graphical representation, or even to create an animation of the 
motion\footnote{Here we use the command \texttt{wxanimate\_draw}, provided by the graphical frontend wxMaxima.}:
\begin{verbatim}

	   surface_hide=true,zrange=[-1.5,1.5],ytics=0.4,xtics=0.4,color=orange,
	   parametric_surface(r*cos(theta),r*sin(theta),subst(t=s,expr),r,0,1,theta,0,2
	),wxanimate_framerate=6$

\end{verbatim}
\begin{center}
\includegraphics[trim=10pt 45pt 10pt 5.5cm, clip,scale=0.625]{wavecilla.pdf}
\end{center}

