Initial commitv0.3

1 files changed, 3679 insertions, 0 deletions

diff --git a/doc/simplesolv-doc.tex b/doc/simplesolv-doc.tex
new file mode 100644
index 0000000..722282a
--- /dev/null
+++ b/doc/simplesolv-doc.tex
@@ -0,0 +1,3679 @@
+%% Copyright 2021 Ian Jauslin
+%% 
+%% Licensed under the Apache License, Version 2.0 (the "License");
+%% you may not use this file except in compliance with the License.
+%% You may obtain a copy of the License at
+%% 
+%%     http://www.apache.org/licenses/LICENSE-2.0
+%% 
+%% Unless required by applicable law or agreed to in writing, software
+%% distributed under the License is distributed on an "AS IS" BASIS,
+%% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+%% See the License for the specific language governing permissions and
+%% limitations under the License.
+
+\documentclass[subsection_in_all]{ian}
+
+\usepackage{code}
+\usepackage{dsfont}
+\usepackage{largearray}
+\usepackage{dlmf}
+
+\def\Eta{H}
+
+\begin{document}
+
+\hbox{}
+\hfil{\bf\Large
+{\tt simplesolv}\par
+\medskip
+\hfil \large v0.3
+}
+
+\vfill
+
+\tableofcontents
+
+\vfill
+\eject
+
+\setcounter{page}1
+\pagestyle{plain}
+
+\indent
+{\tt simplesolv} is a tool to compute the solution of the equations of the ``Simplified approach'' to the repulsive Bose gas introduced in\-~\cite{CJL20,CJL20b,CHe21}.
+This approach provides an approximation to various observables of the ground state of the Hamiltonian
+\begin{equation}
+  H_N=-\frac12\sum_{i=1}^N\Delta_i+\sum_{1\leqslant i\leqslant j\leqslant N}v(|x_i-x_j|)
+\end{equation}
+in three dimensions, with periodic boundary conditions, in the thermodynamic limit $N\to\infty$ at fixed density $\rho$.
+\bigskip
+
+\indent
+{\tt simplesolv} is written in {\tt julia}.
+The source code is located in the {\tt src} directory of this bundle.
+Throughout the documentation, we will refer to the directory containing the {\tt src} directory as the ``installation directory'', and will denote it by the bash variable {\tt\$SIMPLESOLV} (so that the main julia file, for instance, is located at {\tt \$SIMPLESOLV/src/main.jl}).
+\vskip20pt
+
+\section{Basic usage}
+\indent
+Denoting the location of the installation directory by {\tt\$SIMPLESOLV}, {\tt simplesolv} is run by calling
+\begin{code}
+  julia \$SIMPLESOLV/src/main.jl [args] <command>
+\end{code}
+where the optional arguments {\tt [args]} take the form {\tt [-p params] [-U potential] [-M method] [-s savefile]}.
+\bigskip
+
+\indent
+A few commands support multithreaded execution.
+To enable {\tt julia} to run on several processors, it should be run with the {\tt -p} option.
+For example, to run on 8 CPUs, run
+\begin{code}
+  julia -p 8 \$SIMPLESOLV/src/main.jl [args] <command>
+\end{code}
+\bigskip
+
+\indent
+{\tt command} specifies which computation is to be carried out, such as {\tt energy} to compute the ground state energy, or {\tt condensate\_fraction} for the uncondensed fraction.
+The list of available commands depends on the \refname{sec:methods}{{\tt method}} argument, which specifies one of the available methods to solve the equation at hand.
+The available methods are (see section\-~\ref{sec:methods} for further details)
+\begin{itemize}
+  \item \refname{sec:easyeq}{{\tt easyeq}} ({\bf default}) for the Simple or Medium equation, or any interpolation between them, with a soft potential using the Newton algorithm,
+  \item \refname{sec:anyeq}{{\tt anyeq}} for any equation in the ``Simplified approach'' using the Newton algorithm.
+  \item \refname{sec:simpleq-Kv}{{\tt simpleq-Kv}} for the Simple equation using explicit expressions involving $\mathfrak Kv$ (see\-~(\ref{Kv})),
+  \item \refname{sec:simpleq-hardcore}{{\tt simpleq-hardcore}} for the Simple equation with a hard core potential using the Newton algorithm,
+  \item \refname{sec:simpleq-iteration}{{\tt simpleq-iteration}} for the Simple equation with a soft potential using the iteration defined in\-~\cite{CJL20}.
+\end{itemize}
+Each method is described in detail below, along with the list of commands ({\tt command}) and parameters ({\tt params}) compatible with them.
+\bigskip
+
+\indent
+{\tt params} should be a `{\tt;}' separated list of entries, each of which is of the form {\tt key=value}.
+For example {\tt -p "rho=1e-6;v\_a=2"}.
+(Note that you should not end the list of parameters by a `{\tt;}', otherwise {\tt simplesolv} will interpret that as there being an empty parameter entry, which it cannot split into a key and value, and will fail.)
+\bigskip
+
+\indent
+\refname{sec:potentials}{{\tt potential}} specifies which potential $v$ should be used, from the following list (see section\-~\ref{sec:potentials} for further details).
+\begin{itemize}
+  \item \refname{sec:exp}{{\tt exp}} ({\bf default}) for $v(|x|)=ae^{-|x|}$,
+  \item \refname{sec:tent}{{\tt tent}} for $v(|x|)=\mathds 1_{|x|<b}a\frac{2\pi}{3}(1-\frac{|x|}b)^2(\frac{|x|}b+2)$,
+  \item \refname{sec:expcry}{{\tt expcry}} for $v(|x|)=e^{-|x|}-ae^{-b|x|}$,
+  \item \refname{sec:npt}{{\tt npt}} for $v(|x|)=x^2e^{-|x|}$,
+  \item \refname{sec:alg}{{\tt alg}} for $v(|x|)=\frac1{1+\frac14|x|^4}$,
+  \item \refname{sec:algwell}{{\tt algwell}} for $v(|x|)=\frac{1+a|x|^4}{(1+|x|^2)^4}$.
+  \item \refname{sec:exact}{{\tt exact}} for $v(|x|)=\frac{12c( |x|^6b^6(2e-b^2) +b^4|x|^4(9e-7b^2) +4b^2|x|^2(3e-2b^2) +(5e+16b^2))}{(1+b^2|x|^2)^2(4+b^2|x|^2)^2((1+b^2|x|^2)^2-c)}$
+\end{itemize}
+The parameters in the potential can be set using the {\tt params} argument: to set $a$ set {\tt v\_a}, to set $b$ set {\tt v\_b}, to set $c$ set {\tt v\_c}, and to set $e$ set {\tt v\_e}.
+\bigskip
+
+\indent
+{\tt savefile} can be used to accelerate the computation of observables in the \refname{eq:anyeq_compleq}{{\tt compleq}} equation.
+Indeed, as is discussed in section\-~\ref{sec:anyeq}, the computation of \refname{eq:anyeq_compleq}{{\tt compleq}} is based on the computation of a large matrix, which can be pre-computed, saved in a file using the \refname{command:anyeq_save_Abar}{{\tt save\_Abar}} command, and reused by specifying that file in the {\tt savefile} argument.
+
+\section{Methods}\label{sec:methods}
+\indent
+In this section, we describe the different computation methods.
+\bigskip
+
+\subsection{\tt easyeq}\label{sec:easyeq}
+\indent
+This method is used to solve a family of equations, called {\tt easyeq}, that interpolate between the Simple equation and the Medium equation:
+\begin{equation}
+  -\Delta u
+  =
+  v(1-u)-2\rho K+\rho^2 L
+  \label{easyeq}
+\end{equation}
+with
+\begin{equation}
+  K:=\beta_K u\ast S+(1-\beta_K)\frac{2e}\rho u
+  ,\quad
+  L:=\beta_Lu\ast u\ast S+(1-\beta_L)\frac{2e}\rho u\ast u
+  \label{easyeq_KL}
+\end{equation}
+\begin{equation}
+  S:=(1-u)v
+  ,\quad
+  e:=\frac\rho2\int dx\ (1-u(|x|))v(|x|)
+  .
+  \label{easyeq_Se}
+\end{equation}
+for a soft potential $v$ at density $\rho>0$.
+\bigskip
+
+\indent
+\makelink{eq:easyeq_simpleq}{}
+\makelink{eq:easyeq_medeq}{}
+The special choice $\beta_K=\beta_L=0$ is called the Simple equation ({\tt simpleq}), and the choice $\beta_K=\beta_L=1$ is called the Medium equation ({\tt medeq})
+\bigskip
+
+\subsubsection{Usage}
+\indent
+Unless otherwise noted, this method takes the following parameters (specified via the {\tt [-p params]} flag, as a `{\tt;}' separated list of entries).
+\begin{itemize}
+  \item\makelink{param:easyeq_rho}{}
+  {\tt rho} ({\tt Float64}, default: $10^{-6}$): density $\rho$.
+
+  \item\makelink{param:easyeq_tolerance}{}
+  {\tt tolerance} ({\tt Float64}, default: $10^{-11}$): maximal size of final step in Newton iteration.
+
+  \item\makelink{param:easyeq_maxiter}{}
+  {\tt maxiter} ({\tt Int64}, default: 21): maximal number of iterations of the Newton algorithm before giving up.
+
+  \item\makelink{param:easyeq_order}{}
+  {\tt order} ({\tt Int64}, default: 100): order used for all Gauss quadratures (denoted by $N$ below).
+
+  \item\makelink{param:easyeq_minlrho_init}{}
+  {\tt minlrho\_init} ({\tt Float64}, default: {\tt rho}): to initialize the Newton algorithm, we first compute the solution for a smaller $\rho$, {\tt minlrho} is the minimal value for $\log_{10}\rho$ to start this initialization process.
+
+  \item\makelink{param:easyeq_nlrho_init}{}
+  {\tt nlrho\_init} ({\tt Int64}, default: 1): number of steps in the initialization process described above. Set to 1 to disable the incremental initialization process.
+
+  \item\makelink{param:easyeq_bK}{}
+  {\tt bK}, {\tt bL} ({\tt Float64}, default: 1, 1): the values of $\beta_K$ and $\beta_L$.
+
+  \item\makelink{param:easyeq_eq}{}
+  {\tt eq} ({\tt String}, default: ``{\tt simpleq}'', acceptable values: ``{\tt simpleq}'', ``{\tt medeq}''): A shortcut to select either the \refname{eq:easyeq_simpleq}{Simple equation} ($\beta_K=\beta_L=0$) or the \refname{eq:easyeq_medeq}{Medium equation} ($\beta_L=\beta_K=1$). When this option is set, neither {\tt bK} nor {\tt bL} should be set.
+\end{itemize}
+\bigskip
+
+The available {\tt commands} are the following.
+
+\begin{itemize}
+  \item {\tt energy}\makelink{command:easyeq_energy}{}:
+  compute the energy $e$ at a given $\rho$.\par
+  \underline{Output}: [$e$] [Newton error $\epsilon$].
+
+  \item {\tt energy\_rho}\makelink{command:easyeq_energy_rho}{}:
+  compute the energy $e$ as a function of $\rho$.
+  The Newton algorithm is initialized with the hardcore scattering solution (\ref{easyeq_init}) for the lowest $\rho$, and with the previously computed $\rho$ for the larger densities.\par
+  \underline{Disabled parameters}: \refname{param:easyeq_rho}{{\tt rho}}, \refname{param:easyeq_minlrho_init}{{\tt minlrho\_init}} and \refname{param:easyeq_nlrho_init}{{\tt nlrho\_init}}.\par
+  \underline{Extra parameters}:
+  \begin{itemize}
+    \item\makelink{param:easyeq_minlrho}{}
+    {\tt minlrho} ({\tt Float64}, default: $10^{-6}$): minimal value for $\log_{10}\rho$.
+
+    \item\makelink{param:easyeq_maxlrho}{}
+    {\tt maxlrho} ({\tt Float64}, default: $10^{2}$): maximal value for $\log_{10}\rho$.
+
+    \item\makelink{param:easyeq_nlrho}{}
+    {\tt nlrho} ({\tt Int64}, default: $100$): number of values for $\rho$ (spaced logarithmically).
+
+    \item\makelink{param:easyeq_rhos}{}
+    {\tt rhos} ({\tt Array\{Float64\}}, default: $(10^{{\tt minlrho}+\frac{{\tt maxlrho}-{\tt minlrho}}{{\tt nlrho}}n})_n$: list of values for $\rho$, specified as a `{\tt,}' separated list.
+  \end{itemize}
+  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$e$] [Newton error $\epsilon$].
+
+  \item\makelink{command:easyeq_condensate_fraction}{}
+  {\tt condensate\_fraction}: compute the uncondensed fraction $\eta$ at a given $\rho$.\par
+  \underline{Output}: [$\eta$] [Newton error $\epsilon$].
+
+  \item\makelink{command:easyeq_condensate_fraction_rho}{}
+  {\tt condensate\_fraction\_rho}: compute the uncondensed fraction $\eta$ as a function of $\rho$.
+  The Newton algorithm is initialized with the hardcore scattering solution (\ref{easyeq_init}) for the lowest $\rho$, and with the previously computed $\rho$ for the larger densities.\par
+  \underline{Disabled parameters}: same as \refname{command:easyeq_energy_rho}{{\tt energy\_rho}}.\par
+  \underline{Extra parameters}: same as \refname{command:easyeq_energy_rho}{{\tt energy\_rho}}.\par
+  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$\eta$] [Newton error $\epsilon$].
+
+  \item\makelink{command:easyeq_uk}{}
+  {\tt uk}: compute the Fourier transform $\hat u(|k|)$.
+  The values $|k|$ at which $\hat u$ is computed are those coming from the Gauss quadratures, and cannot be set.\par
+  \underline{Output} (one line for each value of $|k|$): [$|k|$] [$\hat u(|k|)$]
+
+  \item\makelink{command:easyeq_ux}{}
+  {\tt ux}: compute $u$ as a function of $|x|$.\par
+  \underline{Extra parameters}:
+  \begin{itemize}
+    \item\makelink{param:easyeq_xmin}{}
+    {\tt xmin} ({\tt Float64}, default: 0): minimum of the range of $|x|$ to be printed.
+
+    \item\makelink{param:easyeq_xmax}{}
+    {\tt xmax} ({\tt Float64}, default: 100): maximum of the range of $|x|$ to be printed.
+
+    \item\makelink{param:easyeq_nx}{}
+    {\tt nx} ({\tt Int64}, default: 100): number of points to print (linearly spaced).
+  \end{itemize}
+  \underline{Output} (one line for each value of $x$): [$|x|$] [$u(|x|)$]
+
+  \item\makelink{command:easyeq_uux}{}
+  {\tt uux}: compute $2u-\rho u\ast u$ as a function of $|x|$.\par
+  \underline{Extra parameters}: Same as \refname{command:easyeq_ux}{{\tt ux}}.\par
+  \underline{Output} (one line for each value of $x$): [$|x|$] [$2u(|x|)-\rho u\ast u(|x|)$]
+
+\end{itemize}
+
+\subsubsection{Description}
+\point{\bf Fourier space formulation.}
+The computation is carried out in Fourier space.
+We take the convention that the Fourier transform of a function $f(|x|)$ is
+\begin{equation}
+  \hat f(|k|)=\int_{\mathbb R^3} dx\ e^{ikx}f(|x|)
+  =\frac{4\pi}{|k|}\int_0^\infty dr\ \frac{\sin(|k|r)}rf(r)
+  .
+\end{equation}
+In Fourier space, (\ref{easyeq}) becomes
+\begin{equation}
+  k^2\hat u=\hat S-2\rho\hat A_K\hat u+\rho^2\hat A_L\hat u^2
+\end{equation}
+\begin{equation}
+  \hat A_K:=\beta_K\hat S+(1-\beta_K)\frac{2e}\rho
+  ,\quad
+  \hat A_L:=\beta_L\hat S+(1-\beta_L)\frac{2e}\rho
+\end{equation}
+with
+\begin{equation}
+  \hat S(|k|)=\int_{\mathbb R^3} dx\ e^{ikx}(1-u(|x|))v(|x|)
+  =\hat v(k)-\hat u\hat\ast\hat v(k)
+  ,\quad
+  \hat f\hat\ast\hat g(|k|):=\int_{\mathbb R^3}\frac{dp}{(2\pi)^3}\ \hat f(|k-p|)\hat g(|p|)
+  .
+\end{equation}
+We write this as a quadratic equation for $\hat u$, and solve it, keeping the solution that decays as $|k|\to\infty$:
+\begin{equation}
+  \rho\hat u(|k|)=
+  \frac{A_K(|k|)}{A_L(|k|)}\left(\xi(|k|)+1-\sqrt{(\xi(|k|)+1)^2-\frac{A_L(|k|)}{A_K^2(|k|)}\hat S(|k|)}\right)
+\end{equation}
+that is,
+\begin{equation}
+  \rho\hat u(|k|)=
+  \frac{\hat S(|k|)}{2A_K(|k|)(\xi(|k|)+1)}\Phi\left({\textstyle\frac{A_L(|k|)}{A_K^2(|k|)}\frac{\hat S(|k|)}{(\xi(|k|)+1)^2}}\right)
+  \label{uk_easyeq}
+\end{equation}
+with
+\begin{equation}
+  \xi(|k|):=\frac{k^2}{2\rho\hat A_K}
+  ,\quad
+  \Phi(x):=\frac{2(1-\sqrt{1-x})}x
+  .
+\end{equation}
+Furthermore, using bipolar coordinates (see lemma\-~\ref{lemma:bipolar}), we write $\hat S$ as
+\begin{equation}
+  \hat S(|k|)=\hat v(|k|)-\frac1{8\pi^3}\int_0^\infty dt\ t\hat u(t)\Eta(|k|,t)
+  ,\quad
+  \Eta(y,t):=\frac{2\pi}y\int_{|y-t|}^{y+t}ds\ s\hat v(s)
+  .
+  \label{easyeq_Eta}
+\end{equation}
+By a simple change of variables,
+\begin{equation}
+  \Eta(y,t)=
+  4\pi\left(\mathds 1_{y>t}\frac{t}y+\mathds 1_{y\leqslant t}\right)
+  \int_0^1 ds\ ((y+t)s+|y-t|(1-s))v((y+t)s+|y-t|(1-s))
+  .
+  \label{Eta}
+\end{equation}
+\bigskip
+
+\point{\bf Evaluating integrals.}
+To compute these integrals numerically, we will use Gauss-Legendre quadratures:
+\begin{equation}
+  \int_0^1 dt\ f(t)\approx\frac12\sum_{i=1}^N w_if\left({\textstyle\frac{r_i+1}2}\right)
+\end{equation}
+where $w_i$ and $r_i$ are the Gauss-Legendre {\it weights} and {\it abcissa}.
+The order $N$ corresponds to the parameter \refname{param:easyeq_order}{{\tt order}}.
+The error made by the quadrature is estimated in appendix\-~\ref{appGL}.
+We compactify the integrals to the interval $(0,1)$: $y:=\frac1{t+1}$:
+\begin{equation}
+  \hat S(|k|)=\hat v(|k|)-\frac1{8\pi^3}\int_0^1 dy\ \frac{(1-y)\hat u(\frac{1-y}y)\Eta(|k|,\frac{1-y}y)}{y^3}
+\end{equation}
+so, using the Gauss-Legendre quadrature, we approximate
+\begin{equation}
+  \hat S(|k|)\approx\hat v(|k|)-\frac1{16\pi^3}\sum_{j=1}^N w_j \frac{(1-y_j)\hat u(\frac{1-y_j}{y_j})\Eta(|k|,\frac{1-y_j}{y_j})}{y_j^3}
+  ,\quad
+  y_j:=\frac{r_j+1}2
+  .
+\end{equation}
+This suggests a natural discretization of Fourier space: let
+\begin{equation}
+  k_i:=\frac{1-r_i}{1+r_i}\equiv\frac{1-y_i}{y_i}
+  .
+  \label{easyeq_ki}
+\end{equation}
+Thus, defining
+\begin{equation}
+  \mathbb U_i:=\rho\hat u(k_i)
+  ,\quad
+  \hat v_i:=\hat v(k_i)
+\end{equation}
+and approximate\-~(\ref{uk_easyeq}):
+\begin{equation}
+  \mathbb U_i=\frac{\mathbb T_i}{2(\mathbb X_i+1)}\Phi\left({\textstyle\mathbb B_i\frac{\mathbb T_i}{(\mathbb X_i+1)^2}}\right)
+  \label{bbU_easyeq}
+\end{equation}
+with
+\begin{equation}
+  \mathbb A_{K,i}:=\beta_K\mathbb S_i+(1-\beta_K)\mathbb E
+  ,\quad
+  \mathbb B_{i}:=\frac{\beta_L\mathbb S_i+(1-\beta_L)\mathbb E}{\mathbb A_{K,i}}
+  ,\quad
+  \mathbb T_i:=\frac{\mathbb S_i}{\mathbb A_{K,i}}
+\end{equation}
+\begin{equation}
+  \mathbb S_i:=\hat v_i-\frac1{16\pi^3\rho}\sum_{j=1}^N w_j\frac{(1-y_j)\mathbb U_j\Eta(k_i,k_j)}{y_j^3}
+  ,\quad
+  \mathbb E:=\hat v(0)-\frac1{16\pi^3\rho}\sum_{j=1}^N w_j\frac{(1-y_j)\mathbb U_j\Eta(0,k_j)}{y_j^3}
+\end{equation}
+\begin{equation}
+  \mathbb X_i:=\frac{k_i^2}{2\rho\mathbb A_{K,i}}
+  .
+\end{equation}
+This is a discrete equation for the vector $(\mathbb U_i)_{i=1}^N$.
+\bigskip
+
+\point{\bf Main algorithm to compute $\mathbb U$.}\par\penalty10000
+\medskip\penalty10000
+\subpoint
+We rewrite\-~(\ref{bbU_easyeq}) as a root finding problem:
+\begin{equation}
+  \Xi_i(\mathbb U):=\mathbb U_i-\frac{\mathbb T_i}{2(\mathbb X_i+1)}\Phi\left({\textstyle B_i\frac{\mathbb T_i}{(\mathbb X_i+1)^2}}\right)
+  =0
+  \label{root_easyeq}
+\end{equation}
+which we solve using the Newton algorithm, that is, we define a sequence of $\mathbb U$'s:
+\begin{equation}
+  \mathbb U^{(n+1)}=\mathbb U^{(n)}-(D\Xi(\mathbb U^{(n)}))^{-1}\Xi(\mathbb U^{(n)})
+\end{equation}
+where $D\Xi$ is the Jacobian of $\Xi$:
+\begin{equation}
+  (D\Xi)_{i,j}:=\frac{\partial\Xi_i}{\partial\mathbb U_j}
+  .
+\end{equation}
+\bigskip
+
+\subpoint
+For small values of $\rho$, we initialize the algorithm with the hardcore scattering solution
+\begin{equation}
+  \hat u_0(k)=\frac{4\pi a_0}{k^2}
+  \label{easyeq_init}
+\end{equation}
+where $a_0$ is the scattering length of the potential $v$ (or an approximation thereof, which need not be very good).
+Thus,
+\begin{equation}
+  \mathbb U^{(0)}_i=\frac{4\pi a_0}{k_i^2}
+  .
+\end{equation}
+This is a good approximation for small $\rho$.
+For larger $\rho$, we choose $\mathbb U^{(0)}$ as the solution of {\tt easyeq} for a slightly smaller $\rho$, and proceed inductively (using the parameters \refname{param:easyeq_minlrho_init}{{\tt minlrho\_init}} and \refname{param:easyeq_nlrho_init}{{\tt nlrho\_init}}).
+\bigskip
+
+\subpoint
+We are left with computing the Jacobian of $\Xi$:
+\begin{equation}
+  \begin{largearray}
+    \frac{\partial\Xi_i}{\partial\mathbb U_j}
+    =\delta_{j,i}
+    -\frac1{2(\mathbb X_i+1)}\left(
+      \partial_j\mathbb T_i-\frac{\mathbb T_i\partial_j\mathbb X_i}{\mathbb X_i+1}
+    \right)\Phi\left({\textstyle\mathbb B_i\frac{\mathbb T_i}{(\mathbb X_i+1)^2}}\right)
+    \\[0.5cm]\hfill
+    -\frac{\mathbb T_i}{2(\mathbb X_i+1)^3}\left(
+      \mathbb B_i\partial_j\mathbb T_i+\mathbb T_i\partial_j\mathbb B_i-2\frac{\mathbb B_i\mathbb T_i\partial_j\mathbb X_i}{\mathbb X_i+1}
+    \right)\partial\Phi\left({\textstyle\mathbb B_i\frac{\mathbb T_i}{(\mathbb X_i+1)^2}}\right)
+  \end{largearray}
+  \label{jacobian_easyeq}
+\end{equation}
+with
+\begin{equation}
+  \partial_j\mathbb B_i=
+  (\beta_L(1-\beta_K)-\beta_K(1-\beta_L))
+  \frac{\mathbb E\partial_j\mathbb S_i-\mathbb S_i\partial_j\mathbb E}{(\beta_K\mathbb S_i+(1-\beta_K)\mathbb E)^2}
+\end{equation}
+\begin{equation}
+  \partial_j\mathbb T_i=
+  (1-\beta_K)\frac{\mathbb E\partial_j\mathbb S_i-\mathbb S_i\partial_j\mathbb E}{(\beta_K\mathbb S_i+(1-\beta_K)\mathbb E)^2}
+  ,\quad
+  \partial_j\mathbb A_{K,i}=\beta_K\partial_j\mathbb S_i+(1-\beta_K)\partial_j\mathbb E
+\end{equation}
+\begin{equation}
+  \partial_j\mathbb S_i:=-\frac1{16\pi^3\rho}w_j\frac{(1-y_j)\Eta(k_i,k_j)}{y_j^3}
+  ,\quad
+  \partial_j\mathbb E:=-\frac1{16\pi^3\rho}\sum_{j=1}^N w_j\frac{(1-y_j)\Eta(0,k_j)}{y_j^3}
+\end{equation}
+\begin{equation}
+  \partial_j\mathbb X_i:=-\frac{k_i^2}{2\rho\mathbb A_{K,i}^2}\partial_j\mathbb A_{K,i}
+  .
+\end{equation}
+\bigskip
+
+\subpoint
+We iterate the Newton algorithm until the Newton relative error $\epsilon$ becomes smaller than the \refname{param:easyeq_tolerance}{{\tt tolerance}} parameter.
+The Newton error is defined as
+\begin{equation}
+  \epsilon=\frac{\|\mathbb U^{(n+1)}-\mathbb U^{(n)}\|_2}{\|\mathbb U^{(n)}\|_2}
+\end{equation}
+where $\|\cdot\|_2$ is the $l_2$ norm.
+The energy thus obtained is
+\begin{equation}
+  e=\frac\rho2\mathds E
+\end{equation}
+the Fourier transform $\hat u$ of the solution is
+\begin{equation}
+  \hat u(k_i)\approx\frac{\mathbb U_i}\rho
+\end{equation}
+and the solution $u$ in real space is obtained by inverting the Fourier transform
+\begin{equation}
+  \begin{array}{>\displaystyle r@{\ }>\displaystyle l}
+    u(|x|)=\int\frac{dk}{8\pi^3}\ e^{-ikx}\hat u(|k|)
+    =&\frac1{2\pi^2|x|}\int_0^1dy\ \frac{(1-y)\sin(|x|\frac{1-y}y)\hat u(\frac{1-y}y)}{y^3}
+    \\[0.5cm]
+    \approx&
+    \frac1{4\pi^2|x|}\sum_{j=1}^Nw_j (1+k_j)^2k_j\sin(|x|k_j)\hat u(k_j)
+    .
+  \end{array}
+  \label{ux}
+\end{equation}
+To compute $2u-\rho u\ast u$, we replace $\hat u$ with $2\hat u-\rho\hat u^2$ in the previous equation.
+\bigskip
+
+\point{\bf Condensate fraction.}
+Finally, to compute the uncondensed fraction, we solve the modified {\tt easyeq} (see\-~\cite{CJL20b})
+\begin{equation}
+  (-\Delta+2\mu)u_\mu=v(1-u_\mu)-2\rho K+\rho^2L
+  \label{easyeq_mu}
+\end{equation}
+where $K$ and $L$ are defined as in\-~(\ref{easyeq_KL})-(\ref{easyeq_Se}) in which $u$ is replaced with $u_\mu$.
+The uncondensed fraction is then
+\begin{equation}
+  \eta=\partial_\mu e|_{\mu=0}
+  =-\frac\rho2\int dx\ v(|x|)\partial_\mu u_\mu(|x|)|_{\mu=0}
+  .
+\end{equation}
+To compute the energy in the presence of the parameter $\mu$, we proceed in the same way as for $\mu=0$, the only difference being that $k^2$ should formally be replaced by $k^2+2\mu$.
+In other words, we consider $\mathbb U_i=u_\mu(|k_i|)$ and define $\Xi(\mathbb U,\mu)$ in the same way as in\-~(\ref{root_easyeq}), except that $\mathbb X_i$ should be replaced by
+\begin{equation}
+  \frac{k_i^2+2\mu}{2\rho\mathbb A_{K,i}}
+  .
+\end{equation}
+We then solve
+\begin{equation}
+  \Xi(\mathbb U,\mu)=0
+  .
+\end{equation}
+By differentiating this identity with respect to $\mu$, we find $\partial_\mu u_\mu$:
+\begin{equation}
+  D\Xi\partial_\mu \mathbb U=-\partial_\mu\Xi
+\end{equation}
+and the uncondensed fraction is
+\begin{equation}
+  \eta
+  =\frac12\int dx\ v(x)(D\Xi)^{-1}\partial_\mu\Xi|_{\mu=0}
+\end{equation}
+thus
+\begin{equation}
+  \eta=\frac1{16\pi^3}\int_0^1 dy\ \frac{(1-y)\Eta(0,\frac{1-y}y)(D\Xi)^{-1}\partial_\mu\Xi|_{\mu=0}(\frac{1-y}y)}{y^3}
+\end{equation}
+which we approximate using a Gauss-Legendre quadrature:
+\begin{equation}
+  \eta\approx\frac1{32\pi^3}\sum_{j=1}^N w_j(1-k_j)^2k_j\Eta(0,k_j)\sum_{i=1}^N(D\Xi)^{-1}_{j,i}\partial_\mu\Xi_i|_{\mu=0}
+  .
+  \label{eta_easyeq}
+\end{equation}
+We then compute, using\-~(\ref{jacobian_easyeq}),
+\begin{equation}
+  \partial_\mu\Xi_i|_{\mu=0}
+  =
+  \frac{\mathbb T_i}{2\rho\mathbb A_{K,i}(\mathbb X_i+1)^2}\Phi\left({\textstyle\mathbb B_i\frac{\mathbb T_i}{(\mathbb X_i+1)^2}}\right)
+  +\frac{\mathbb B_i\mathbb T_i^2}{\rho\mathbb A_{K,i}(\mathbb X_i+1)^4}\partial\Phi\left({\textstyle\mathbb B_i\frac{\mathbb T_i}{(\mathbb X_i+1)^2}}\right)
+  .
+\end{equation}
+
+\subsection{\tt anyeq}\label{sec:anyeq}
+\indent
+This method is used to solve any of the equations in the Simplified approach.
+Specifically, it solves the equation
+\begin{equation}
+  -\Delta u
+  =
+  v(1-u)-2\rho K+\rho^2 L
+  \label{anyeq}
+\end{equation}
+with
+\begin{equation}
+  K=\gamma_K(1-\alpha_Ku)\left(\beta_K u\ast S+(1-\beta_K)\frac{2e}\rho u\right)
+  \label{anyeq_K}
+\end{equation}
+and
+\begin{equation}
+  L:=L_1+L_2+L_3
+\end{equation}
+\begin{equation}
+  L_1:=(1-\alpha_{L,1}u)\left(\beta_{L,1}u\ast u\ast S+(1-\beta_{L,1})\frac{2e}\rho u\ast u\right)
+\end{equation}
+\begin{equation}
+  L_2:=-\gamma_{L,2}(1-\alpha_{L,2}u)\left(\beta_{L,2}2u\ast(u(u\ast S))+(1-\beta_{L,2})\frac{4e}\rho u\ast u^2\right)
+\end{equation}
+\begin{equation}
+  L_3:=\gamma_{L,3}(1-\alpha_{L,3}u)\left(\beta_{L,3}\frac12\int dydz\ u(y)u(z-x)u(z)u(y-x)S(z-y)+(1-\beta_{L,3})\frac e\rho u^2\ast u^2\right)
+\end{equation}
+\begin{equation}
+  e=\frac\rho2\int dx\ (1-u(|x|))v(|x|)
+  .
+  \label{anyeq_e}
+\end{equation}
+The parameters $\alpha_\cdot$, $\beta_\cdot$ and $\gamma_\cdot$ can be set to turn\-~(\ref{anyeq}) into any of the approximations of the Simplified approach.
+For ease of use, there are several predefined equations, given in the following table.
+\bigskip
+
+\def\arraystretch{1.5}
+\setlength\tabcolsep{5pt}
+\makelink{table:anyeq_eqs}{}
+\hfil\begin{tabular}{|r|c|c|c|c|c|c|c|c|c|c|c|c|}
+  \hline
+  &$\alpha_K$ &$\beta_K$ &$\gamma_K$ &$\alpha_{L,1}$ &$\beta_{L,1}$ &$\alpha_{L,2}$ &$\beta_{L,2}$ &$\gamma_{L,2}$ &$\alpha_{L,3}$ &$\beta_{L,3}$ &$\gamma_{L,3}$\\
+  \hline
+  \makelink{eq:anyeq_compleq}{}{\tt compleq}&1&1&1&1&1&1&1&1&1&1&1\\
+  \makelink{eq:anyeq_bigeq}{}  {\tt bigeq  }&1&1&1&1&1&1&1&1&-&-&0\\
+  \makelink{eq:anyeq_fulleq}{} {\tt fulleq }&1&1&1&1&1&1&1&1&1&0&1\\
+  \makelink{eq:anyeq_medeq}{}  {\tt medeq  }&0&1&1&0&1&-&-&0&-&-&0\\
+  \makelink{eq:anyeq_simpleq}{}{\tt simpleq}&0&0&1&0&0&-&-&0&-&-&0\\
+  \hline
+\end{tabular}
+\bigskip
+
+Note that there is no $\gamma_{L,1}$, whose computation would be rather different.
+Note, in addition, that {\tt simpleq} and {\tt medeq} coincide with their definitions in\-~(\ref{easyeq}).
+The method used to solve this equation is very different from \refname{sec:easyeq}{{\tt easyeq}}, and is significantly longer to run.
+\bigskip
+
+\subsubsection{Usage}
+\indent
+Unless otherwise noted, this method takes the following parameters (specified via the {\tt [-p params]} flag, as a `{\tt;}' separated list of entries).
+\begin{itemize}
+  \item\makelink{param:anyeq_rho}{}
+  {\tt rho} ({\tt Float64}, default: $10^{-6}$): density $\rho$.
+
+  \item\makelink{param:anyeq_tolerance}{}
+  {\tt tolerance} ({\tt Float64}, default: $10^{-11}$): maximal size of final step in Newton iteration.
+
+  \item\makelink{param:anyeq_maxiter}{}
+  {\tt maxiter} ({\tt Int64}, default: 21): maximal number of iterations of the Newton algorithm before giving up.
+
+  \item\makelink{param:anyeq_P}{}
+  {\tt P} ({\tt Int64}, default: 11): order of all Chebyshev polynomial expansions (denoted by $P$ below).
+
+  \item\makelink{param:anyeq_N}{}
+  {\tt N} ({\tt Int64}, default: 12): order of all Gauss quadratures (denoted by $N$ below).
+
+  \item\makelink{param:anyeq_J}{}
+  {\tt J} ({\tt Int64}, default: 10): number of splines (denoted by $J$ below).
+
+  \item\makelink{param:anyeq_minlrho_init}{}
+  {\tt minlrho\_init} ({\tt Float64}, default: \refname{param:anyeq_rho}{{\tt rho}}): we initialize the Newton algorithm using the solution of \refname{eq:easyeq_medeq}{{\tt medeq}}, computed using the methods in \refname{sec:easyeq}{{\tt easyeq}}. This option is passed to the underlying \refname{sec:easyeq}{{\tt easyeq}} routine.
+
+  \item\makelink{param:anyeq_nlrho_init}{}
+  {\tt nlrho\_init} ({\tt Int64}, default: 1): this option is passed to the underlying \refname{sec:easyeq}{{\tt easyeq}} routine to initialize the Newton algorithm.
+
+  \item\makelink{param:anyeq_aK}{}
+  {\tt aK}, {\tt bK}, {\tt gK}, {\tt aL1}, {\tt bL1}, {\tt aL2}, {\tt bL2}, {\tt gL2}, {\tt aL3}, {\tt bL3}, {\tt gL3} ({\tt Float64}, default: 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0): the values of $\alpha_K$, $\beta_K$, $\gamma_K$, $\alpha_{L,1}$, $\beta_{L,1}$, $\alpha_{L,2}$, $\beta_{L,2}$, $\gamma_{L,2}$, $\alpha_{L,3}$, $\beta_{L,3}$, $\gamma_{L,3}$.
+
+  \item\makelink{param:anyeq_eq}{}
+  {\tt eq} ({\tt String}, default: ``{\tt bigeq}'', acceptable values: ``{\tt compleq}'', ``{\tt bigeq}'', ``{\tt fulleq}'', ``{\tt medeq}'', ``{\tt simpleq}''): A shortcut to select any of the equations defined in the \refname{table:anyeq_eqs}{table above}. When this option is set, none of
+  \refname{param:anyeq_aK}{\tt aK}, \refname{param:anyeq_aK}{\tt bK}, \refname{param:anyeq_aK}{\tt gK}, \refname{param:anyeq_aK}{\tt aL1}, \refname{param:anyeq_aK}{\tt bL1}, \refname{param:anyeq_aK}{\tt aL2}, \refname{param:anyeq_aK}{\tt bL2}, \refname{param:anyeq_aK}{\tt gL2}, \refname{param:anyeq_aK}{\tt aL3}, \refname{param:anyeq_aK}{\tt bL3}, \refname{param:anyeq_aK}{\tt gL3} should be set.
+\end{itemize}
+\bigskip
+
+The available {\tt commands} are the following.
+
+\begin{itemize}
+  \item\makelink{command:anyeq_energy}{}
+  {\tt energy}: compute the energy $e$ at a given $\rho$.\par
+  \underline{Output}: [$e$] [Newton error $\epsilon$].
+
+  \item\makelink{command:anyeq_energy_rho}{}
+  {\tt energy\_rho}: compute the energy $e$ as a function of $\rho$.
+  The Newton algorithm is initialized with the solution of \refname{eq:easyeq_medeq}{{\tt medeq}}.\par
+  \underline{Disabled parameters}: \refname{param:anyeq_rho}{{\tt rho}}, \refname{param:anyeq_minlrho_init}{{\tt minlrho\_init}} and \refname{param:anyeq_nlrho_init}{{\tt nlrho\_init}}.\par
+  \underline{Extra parameters}:
+  \begin{itemize}
+    \item\makelink{param:anyeq_minlrho}{}
+    {\tt minlrho} ({\tt Float64}, default: $10^{-6}$): minimal value for $\log_{10}\rho$.
+
+    \item\makelink{param:anyeq_maxlrho}{}
+    {\tt maxlrho} ({\tt Float64}, default: $10^{2}$): maximal value for $\log_{10}\rho$.
+
+    \item\makelink{param:anyeq_nlrho}{}
+    {\tt nlrho} ({\tt Int64}, default: $100$): number of values for $\rho$ (spaced logarithmically).
+
+    \item\makelink{param:anyeq_rhos}{}
+    {\tt rhos} ({\tt Array\{Float64\}}, default: $(10^{{\tt minlrho}+\frac{{\tt maxlrho}-{\tt minlrho}}{{\tt nlrho}}n})_n$: list of values for $\rho$, specified as a `{\tt,}' separated list.
+  \end{itemize}
+  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$e$] [Newton error $\epsilon$].\par
+  \underline{Multithread support}: yes, different values of $\rho$ split up among workers.
+
+  \item\makelink{command:anyeq_energy_rho_init_prevrho}{}
+  {\tt energy\_rho\_init\_prevrho}: compute the energy $e$ as a function of $\rho$.
+  The Newton algorithm is initialized with the solution of \refname{eq:easyeq_medeq}{{\tt medeq}} for the lowest $\rho$, and with the previously computed $\rho$ for the larger densities.\par
+  \underline{Disabled parameters}: same as \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
+  \underline{Extra parameters}: same as \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
+  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$e$] [Newton error $\epsilon$].
+
+  \item\makelink{command:anyeq_energy_rho_init_nextrho}{}
+  {\tt energy\_rho\_init\_nextrho}: same as \refname{command:anyeq_energy_rho_init_prevrho}{{\tt energy\_rho\_init\_prevrho}} except that the energy is computed for decreasing densities instead of increasing ones.
+  The Newton algorithm is initialized with the solution of \refname{eq:easyeq_medeq}{{\tt medeq}} for the largest $\rho$, and with the previously computed $\rho$ for the smaller densities.\par
+  \underline{Disabled parameters}: same as \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
+  \underline{Extra parameters}: same as \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
+  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$e$] [Newton error $\epsilon$].
+
+  \item\makelink{command:anyeq_condensate_fraction}{}
+  {\tt condensate\_fraction}: compute the uncondensed fraction $\eta$ at a given $\rho$.\par
+  \underline{Output}: [$\eta$] [Newton error $\epsilon$].
+
+  \item\makelink{command:anyeq_condensate_fraction_rho}{}
+  {\tt condensate\_fraction\_rho}: compute the uncondensed fraction $\eta$ as a function of $\rho$.
+  The Newton algorithm is initialized with the solution of \refname{eq:easyeq_medeq}{{\tt medeq}}.\par
+  \underline{Disabled parameters}: same as \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
+  \underline{Extra parameters}: same as \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
+  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$\eta$] [Newton error $\epsilon$].\par
+  \underline{Multithread support}: yes, different values of $\rho$ split up among workers.
+
+  \item\makelink{command:anyeq_uk}{}
+  {\tt uk}: compute the Fourier transform $\hat u(|k|)$.
+  The values $|k|$ at which $\hat u$ is computed are those coming from the Gauss quadratures, and cannot be set.\par
+  \underline{Output} (one line for each value of $|k|$): [$|k|$] [$\hat u(|k|)$]
+
+  \item\makelink{command:anyeq_ux}{}
+  {\tt ux}: compute $u$ as a function of $|x|$.\par
+  \underline{Extra parameters}:
+  \begin{itemize}
+    \item\makelink{param:anyeq_xmin}{}
+    {\tt xmin} ({\tt Float64}, default: 0): minimum of the range of $|x|$ to be printed.
+
+    \item\makelink{param:anyeq_xmax}{}
+    {\tt xmax} ({\tt Float64}, default: 100): maximum of the range of $|x|$ to be printed.
+
+    \item\makelink{param:anyeq_nx}{}
+    {\tt nx} ({\tt Int64}, default: 100): number of points to print (linearly spaced).
+
+  \end{itemize}
+  \underline{Output} (one line for each value of $x$): [$|x|$] [$\mathcal Re(u(|x|))$] [$\mathcal Im(u(|x|))$]
+
+  \item\makelink{command:anyeq_momentum_distribution}{}
+  {\tt momentum\_distribution}: compute the momentum distribution $\mathfrak M(|k|)$ at a given $\rho$.\par
+  \underline{Output} (one line for each value of $|k|$): [$|k|$] [$\mathfrak M(|k|)$]
+
+  \item\makelink{command_anyeq_2pt}{}
+  {\tt 2pt}: compute the two-point correlation function $C_2(|x|)$ at a given $\rho$.\par
+  \underline{Extra parameters}: same as \refname{command:anyeq_ux}{{\tt ux}}, plus\par
+  \begin{itemize}
+    \item\makelink{param:anyeq_window_L}{}
+    {\tt window\_L} ({\tt Float64}, default: $10^3$): size of the Hann window used to numerically invert the Fourier transform in the computation of the tow-point correlation function, see\-~(\ref{hann}).
+  \end{itemize}
+  \underline{Output} (one line for each value of $|x|$): [$|x|$] [$C_2(|x|)$]
+
+  \item \makelink{command:anyeq_save_Abar}{}
+  {\tt save\_Abar}: compute the matrix $\bar A$. This matrix is used to compute observables for \refname{eq:anyeq_compleq}{{\tt compleq}}. This command is useful to output the value of $\bar A$ to a file once and for all, and use this file to run commands without recomputing $\bar A$.\par
+  \underline{Disabled parameters}: \refname{param:anyeq_rho}{{\tt rho}}, \refname{param:anyeq_tolerance}{{\tt tolerance}}, \refname{param:anyeq_maxiter}{{\tt maxiter}}, \refname{param:anyeq_minlrho_init}{{\tt minlrho\_init}}, \refname{param:anyeq_nlrho}{{\tt nlrho}}.\par
+  \underline{Output}: [$\bar A$] (the output is not designed to be human-readable; it is obtained through nested {\tt for} loops; for details, see the code).\par
+  \underline{Multithread support}: yes, the first indices are split up among workers, which produces $NJ$ jobs.
+
+\end{itemize}
+
+\subsubsection{Description}
+\point{\bf Fourier space formulation.}
+The computation is carried out in Fourier space.
+We take the convention that the Fourier transform of a function $f(|x|)$ is
+\begin{equation}
+  \hat f(|k|)=\int_{\mathbb R^3} dx\ e^{ikx}f(|x|)
+  =\frac{4\pi}{|k|}\int_0^\infty dr\ \frac{\sin(|k|r)}rf(r)
+  .
+  \label{anyeq_fourier}
+\end{equation}
+We define a Fourier-space convolution:
+\begin{equation}
+  \hat f\hat\ast\hat g(|k|):=\int_{\mathbb R^3}\frac{dp}{(2\pi)^3}\ \hat f(|k-p|)\hat g(|p|)
+  .
+  \label{anyeq_hatast}
+\end{equation}
+In Fourier space, (\ref{anyeq}) becomes
+\begin{equation}
+  k^2\hat u(|k|)=\hat S(|k|)-2\rho\hat K(|k|)+\rho^2\hat L(|k|)
+\end{equation}
+with
+\begin{equation}
+  \hat S(|k|)=\int_{\mathbb R^3} dx\ e^{ikx}(1-u(|x|))v(|x|)
+  =\hat v(k)-\hat u\hat\ast\hat v(k)
+\end{equation}
+\begin{equation}
+  \rho\hat K=
+  \gamma_K(\beta_K\rho\hat S+(1-\beta_K)2e)\hat u
+  -
+  \gamma_K\alpha_K\hat u\hat\ast((\beta_K\rho\hat S+(1-\beta_K)2e)\hat u)
+\end{equation}
+\begin{equation}
+  \rho\hat L_1=
+  (\beta_{L,1}\rho\hat S+(1-\beta_{L,1})2e)\hat u^2
+  -
+  \alpha_{L,1}\hat u\hat\ast((\beta_{L,1}\rho\hat S+(1-\beta_{L,1})2e)\hat u^2)
+\end{equation}
+\begin{equation}
+  \rho\hat L_2=
+  -\gamma_{L,2}(\beta_{L,2}2\rho(\hat u\hat\ast(\hat u\hat S))+(1-\beta_{L,2})4e\hat u\hat\ast\hat u)\hat u
+  +
+  \gamma_{L,2}\alpha_{L,2}\hat u\hat\ast((\beta_{L,2}2\rho(\hat u\hat\ast(\hat u\hat S))+(1-\beta_{L,2})4e\hat u\hat\ast\hat u)\hat u)
+\end{equation}
+\begin{equation}
+  \rho\hat L_3=
+  \gamma_{L,3}(1-\beta_{L,3})e(\hat u\hat\ast\hat u)^2
+  -
+  \gamma_{L,3}\alpha_{L,3}(1-\beta_{L,3})\hat u\hat\ast(e(\hat u\hat\ast\hat u)^2)
+  +
+  \gamma_{L,3}\beta_{L,3}\hat l_3
+  -
+  \frac1{2\rho}\gamma_{L,3}\alpha_{L,3}\beta_{L,3}\hat u\hat\ast\hat l_3
+\end{equation}
+with
+\begin{equation}
+  \hat l_3(|k|):=
+  \frac1{\rho^2}\int\frac{dq}{(2\pi)^3}\frac{dq'}{(2\pi)^3}\ \rho\hat u(|k-q'|)\rho\hat u(|q|)\rho\hat u(|q'|)\rho\hat u(|k-q|)\hat S(|q'-q|)
+  \label{l3}
+\end{equation}
+Therefore,
+\begin{equation}
+  \rho^2\hat u(|k|)^2\sigma_{L,1}(|k|)
+  -\left(\frac{k^2}\rho+2\sigma_K(|k|)+2f_1(|k|)\right)\rho\hat u(|k|)
+  +\left(\hat S(|k|)+\frac12f_2(|k|)+G(|k|)\right)
+  =
+  0
+\end{equation}
+with
+\begin{equation}
+  \sigma_K:=\frac1{\rho}\gamma_K(\beta_K\rho\hat S+(1-\beta_K)2e)
+  ,\quad
+  \sigma_{L,1}:=\frac1{\rho}(\beta_{L,1}\rho\hat S+(1-\beta_{L,1})2e)
+\end{equation}
+\begin{equation}
+  f_1:=\frac1{\rho^2}\gamma_{L,2}(\beta_{L,2}\rho(\rho\hat u\hat\ast(\rho\hat u\hat S))+(1-\beta_{L,2})2e\rho\hat u\hat\ast\rho\hat u)
+\end{equation}
+\begin{equation}
+  f_2:=\gamma_{L,3}(1-\beta_{L,3})\frac{2e}{\rho^3}(\rho\hat u\hat\ast\rho\hat u)^2
+  +\gamma_{L,3}\beta_{L,3}\hat l_3
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    G:=
+    \frac2{\rho^2}\gamma_K\alpha_K\rho\hat u\hat\ast((\beta_K\rho\hat S+(1-\beta_K)2e)\rho\hat u)
+    -
+    \frac1{\rho^2}\alpha_{L,1}\rho\hat u\hat\ast((\beta_{L,1}\rho\hat S+(1-\beta_{L,1})2e)\rho\hat u^2)
+    \\[0.5cm]\hfill
+    +
+    \frac2{\rho}\alpha_{L,2}\rho\hat u\hat\ast(f_1\rho\hat u)
+    -
+    \frac1{2\rho}\alpha_{L,3}\rho\hat u\hat\ast f_2
+    .
+  \end{largearray}
+\end{equation}
+Therefore,
+\begin{equation}
+  \rho\hat u=1+\xi
+  -\sqrt{(1+\xi)^2-(1+\zeta)}
+  \label{u}
+\end{equation}
+with
+\begin{equation}
+  \xi:=\frac1{\sigma_{L,1}}\left(\sigma_K-\sigma_{L,1}+\frac{k^2}{2\rho}+f_1\right)
+  ,\quad
+  \zeta:=\frac1{\sigma_{L,1}}\left(\hat S-\sigma_{L,1}+\frac12f_2+G\right)
+\end{equation}
+We rewrite\-~(\ref{u}) as
+\begin{equation}
+  U=
+  \frac{1+\zeta}{2(\xi+1)}\ \Phi\left({\textstyle\frac{1+\zeta}{(\xi+1)^2}}\right)
+  \label{anyeq_hatu}
+\end{equation}
+with
+\begin{equation}
+  U:=\rho\hat u
+\end{equation}
+\begin{equation}
+  \Phi(x):=\frac{2(1-\sqrt{1-x})}x
+  .
+\end{equation}
+\bigskip
+
+\point{\bf Evaluating integrals.}
+To evaluate integrals numerically, we will split integration intervals over splines and use Gauss-Legendre quadratures.
+More specifically, to compute integrals of the form
+\begin{equation}
+  \int \frac{dk}{(2\pi)^3}\ f(U(k),k)
+  =
+  \frac1{2\pi^2}\int_0^\infty dt\ t^2f(U(t),t)
+\end{equation}
+we first compactify space to $(-1,1]$ by changing variables to $\tau=\frac{1-t}{1+t}$:
+\begin{equation}
+  \int \frac{dk}{(2\pi)^3}\ f(U(k),k)
+  =
+  \frac1{\pi^2}\int_{-1}^1 d\tau\ \frac{(1-\tau)^2}{(1+\tau)^4}f({\textstyle U(\frac{1-\tau}{1+\tau}),\frac{1-\tau}{1+\tau}})
+  .
+\end{equation}
+We then split $(-1,1]$ into $J$ sub-intervals (given by the parameter \refname{param:anyeq_J}{{\tt J}}), called {\it splines}:
+\begin{equation}
+  (-1,1]=\bigcup_{j=0}^{J-1}(\tau_j,\tau_{j+1}]
+  .
+\end{equation}
+The $\tau$ are taken to be equally spaced, but the code is designed in such a way that this could be changed easily in the future:
+\begin{equation}
+  \tau_j=-1+\frac{2j}J
+  .
+\end{equation}
+In these terms,
+\begin{equation}
+  \int \frac{dk}{(2\pi)^3}\ f(U(k),k)
+  =
+  \sum_{l=0}^{J-1}
+  \int_{\tau_l}^{\tau_{l+1}} d\tau\ \frac{(1-\tau)^2}{(1+\tau^4}g(U({\textstyle\frac{1-\tau}{1+\tau},\frac{1-\tau}{1+\tau}})
+  .
+\end{equation}
+We then change variables to $r$:
+\begin{equation}
+  \tau=
+  -\frac{\tau_{l+1}-\tau_l}2\sin({\textstyle\frac{\pi r}2})+\frac{\tau_{l+1}+\tau_l}2
+  =:\frac{1-\vartheta_l(r)}{1+\vartheta_l(r)}
+  \label{chebyshevvars}
+\end{equation}
+(the reason for this specific change of variables will become clear at the end of this paragraph and in the next one)
+and find
+\begin{equation}
+  \int \frac{dk}{(2\pi)^3}\ f(U(k),k)
+  =
+  \sum_{l=0}^{J-1}
+  \frac{\tau_{l+1}-\tau_l}{16\pi}\int_{-1}^1 dr\ \cos({\textstyle\frac{\pi r}2})(1+\vartheta_l(r))^2\vartheta_l^2(r)f(U(\vartheta_l(r)),\vartheta_l(r))
+  .
+\end{equation}
+We then approximate the integral using a Gauss-Legendre quadrature of order $N$ (see appendix\-~\ref{appGL}):
+\begin{equation}
+  \int \frac{dk}{(2\pi)^3}\ f(k)
+  \approx
+  \sum_{l=0}^{J-1}
+  \frac{\tau_{l+1}-\tau_l}{16\pi}
+  \sum_{j=1}^Nw_j
+  \cos({\textstyle\frac{\pi x_j}2})(1+\vartheta_l(x_j))^2\vartheta_l^2(x_j)f(\mathbb U_{l,j},\vartheta_l(x_j))
+  \label{anyeq_integration}
+\end{equation}
+with
+\begin{equation}
+  \mathbb U_{l,j}:=U(k_{l,j})
+  ,\quad
+  k_{l,j}:=\frac{2+(\tau_{l+1}-\tau_l)\sin(\frac{\pi x_j}2)-(\tau_{l+1}+\tau_\alpha)}{2-(\tau_{l+1}-\tau_l)\sin(\frac{\pi x_j}2)+(\tau_{l+1}+\tau_l)}
+  .
+\end{equation}
+The choice of the change of variables\-~(\ref{chebyshevvars}) is made so that $U$ is evaluated at $k_{l,j}$, which are the momenta that appear naturally in the Chebyshev polynomial expansion below\-~(\ref{klj}).
+In this way, we can compute arbitrary integrals of functions of $U$ by just computing the values of $\mathbb U_{l,j}$.
+\bigskip
+
+\point{\bf Chebyshev polynomial expansion.}
+In the case of \refname{sec:easyeq}{{\tt easyeq}}, we saw that using Gauss quadratures reduced the computation to evaluating $U$ at a fixed set of discrete momenta.
+This is not the case here, due to the presence of more complicated convolutions of $U$.
+Instead, we will approximate $U$ using polynomial approximations.
+We will now discuss this approximation for an arbitrary function $a(|k|)$, as we will need it for other functions than simply $U$.
+\bigskip
+
+\subpoint
+As in the previous paragraph, we compactify $[0,\infty)$ to $(-1,1]$ via the change of variables $r\mapsto\frac{1-r}{1+r}$, and split the interval into splines.
+Inside each spline, we use a Chebyshev polynomial expansion.
+However, we need to be careful in the first spline, which encapsulates the behavior at $|k|\to\infty$: $U$ decays to 0 as $|k|\to\infty$, and this behavior cannot be reproduced by a polynomial expansion.
+To avoid this, if $a\sim|k|^{-\nu_a}$, we expand $|k|^{\nu_a}a$ instead of $a$.
+Actually, to simplify expressions in the presence of the compactification, we will expand $2^{-\nu_a}(1+|k|)^{\nu_a}a$, which makes no conceptual difference.
+In the case of $a=U$, we will use $\nu=2$, which we know holds for the Simple equation\-~\cite{CJL20}.
+Putting all this together, we write, for $\tau\in(-1,1]$, if $|k|\equiv\frac{1-\tau}{1+\tau}$,
+\begin{equation}
+  2^{-\nu_a}(1+|k|)^{\nu_a}a(|k|)\equiv
+  \frac1{(1+\tau)^{\nu_a}}a({\textstyle\frac{1-\tau}{1+\tau}})
+  =\sum_{l=0}^{J-1}\mathds 1_{\tau_l<\tau\leqslant\tau_{l+1}}\sum_{n=0}^\infty \mathbf F_{l,n}^{(\nu_a)}(a) T_n({\textstyle\frac{2\tau-(\tau_l+\tau_{l+1})}{\tau_{l+1}-\tau_l}})
+  \label{a_chebyshev_spline}
+\end{equation}
+in which $\mathds 1_{\tau_l<\tau\leqslant\tau_{l+1}}\in\{0,1\}$ is equal to 1 if and only if $\tau_l<\tau\leqslant\tau_{l+1}$, and
+\begin{equation}
+  \mathbf F_{l,n}^{(\nu_a)}(a):=
+  \frac{2-\delta_{n,0}}{\pi}
+  \int_0^{\pi} d\theta\ 
+  \frac{\cos(n\theta)}{(1+\frac{\tau_{l+1}-\tau_l}2\cos(\theta)+\frac{\tau_{l+1}+\tau_l}2)^{\nu_a}}
+  a({\textstyle\frac{2-(\tau_{l+1}-\tau_l)\cos(\theta)-(\tau_{l+1}+\tau_\alpha)}{2+(\tau_{l+1}-\tau_l)\cos(\theta)+(\tau_{l+1}+\tau_l)}})
+  \label{bfF}
+\end{equation}
+(see appendix\-~\ref{appChebyshev} for a discussion of the Chebyshev polynomial expansion and the error of its truncations).
+\bigskip
+
+\subpoint
+In order to compute an approximation for $a$ using\-~(\ref{a_chebyshev_spline}), we will truncate the sum over $n$ to a finite value $P$ (given by the parameter \refname{param:anyeq_P}{{\tt P}}).
+In addition, to compute the integral in\-~(\ref{bfF}), we will use a Gauss-Legendre quadrature of order $N$ (given by the parameter \refname{param:anyeq_N}{{\tt N}}), see appendix\-~\ref{appGL}:
+\begin{equation}
+  \mathbf F_{l,n}^{(\nu_a)}(a)\approx\mathfrak F_{l,n}^{(\nu_a)}(\mathfrak a):=
+  \frac{2-\delta_{n,0}}2\sum_{j=1}^N w_j \cos({\textstyle\frac{n\pi(1+x_j)}2})\frac{\mathfrak a_{l,j}}{(1-\frac{\tau_{l+1}-\tau_l}2\sin(\frac{\pi x_j}2)+\frac{\tau_{l+1}+\tau_l}2)^{\nu_a}}
+  \label{frakF}
+\end{equation}
+with
+\begin{equation}
+  \mathfrak a_{l,j}:=a(k_{l,j})
+  ,\quad
+  k_{l,j}:=\frac{2+(\tau_{l+1}-\tau_l)\sin(\frac{\pi x_j}2)-(\tau_{l+1}+\tau_\alpha)}{2-(\tau_{l+1}-\tau_l)\sin(\frac{\pi x_j}2)+(\tau_{l+1}+\tau_l)}
+  \label{klj}
+\end{equation}
+and $(x_j,w_j)$ are the abcissa and weights for Gauss-Legendre quadratures.
+\bigskip
+
+\subpoint
+All in all, we approximate
+\begin{equation}
+  a({\textstyle\frac{1-\tau}{1+\tau}})
+  \approx(1+\tau)^{\nu_a}\sum_{l=0}^{J-1}\mathds 1_{\tau_l<\tau\leqslant\tau_{l+1}}\sum_{n=0}^P \mathfrak F_{l,n}^{(\nu_a)}(\mathfrak a) T_n({\textstyle\frac{2\tau-(\tau_l+\tau_{l+1})}{\tau_{l+1}-\tau_l}})
+  \label{approxchebyshev}
+\end{equation}
+with $\mathfrak F$ defined in\-~(\ref{frakF}).
+Furthermore, using the Chebyshev polynomial expansion and Gauss-Legendre quadratures, we can compute all the observables we are interested in by computing $\mathbb U_{l,j}\equiv U(k_{l,j})$.
+With this in mind, we will represent the function $U$ as a vector of dimension $NJ$ whose components are $\mathbb U_{l,j}$.
+\bigskip
+
+\point{\bf Convolutions.}
+Using the Chebyshev polynomial approximation, we can compute convolutions as follows.
+\medskip
+
+\subpoint
+First, we rewrite the convolution as a two-dimensional integral, using bipolar coordinates (see lemma\-~\ref{lemma:bipolar}):
+\begin{equation}
+  (a\hat\ast b)(|k|)=\frac1{4\pi^2|k|}\int_0^\infty dt\ ta(t)\int_{||k|-t|}^{|k|+t}ds\ sb(s)
+\end{equation}
+We change variables to compactify the integrals $\tau=\frac{1-t}{1+t}$, $\sigma=\frac{1-s}{1+s}$:
+\begin{equation}
+  (a\hat\ast b)(|k|)=\frac1{4\pi^2|k|}\int_{-1}^1 d\tau\ \frac{2(1-\tau)}{(1+\tau)^3}a({\textstyle\frac{1-\tau}{1+\tau}})\int_{\alpha_-(|k|,\tau)}^{\alpha_+(|k|,\tau)}d\sigma\ \frac{2(1-\sigma)}{(1+\sigma)^3}b({\textstyle\frac{1-\sigma}{1+\sigma}})
+  \label{convvar}
+\end{equation}
+with
+\begin{equation}
+  \alpha_-(|k|,\tau):=\frac{1-|k|-\frac{1-\tau}{1+\tau}}{1+|k|+\frac{1-\tau}{1+\tau}}
+  ,\quad
+  \alpha_+(|k|,\tau):=\frac{1-||k|-\frac{1-\tau}{1+\tau}|}{1+||k|-\frac{1-\tau}{1+\tau}|}
+  .
+  \label{anyeq_alpha}
+\end{equation}
+Therefore, using the approximation\-~(\ref{approxchebyshev}), if $\mathfrak a_{l,j}:=a(k_{l,j})$ and $\mathfrak b_{l,j}:=b(k_{l,j})$ (\ref{klj}),
+\begin{equation}
+  (a\hat\ast b)(|k_{l,j}|)\approx
+  (\mathfrak a\odot \mathfrak b)_{l,j}:=
+  \frac1{4\pi^2|k_{l,j}|}
+  \sum_{n,m=0}^P
+  \sum_{l,l'=0}^{J-1}
+  \mathfrak F_{l,n}^{(\nu_a)}(\mathfrak a)\mathfrak F_{l',m}^{(\nu_b)}(\mathfrak b)
+  A_{l,n;l',m}^{(\nu_a,\nu_b)}(|k_{l,j}|)
+  \label{odot}
+\end{equation}
+with
+\begin{equation}
+  \begin{largearray}
+    A_{l,n;l',m}^{(\nu_a,\nu_b)}(|k|):=
+    \int_{\tau_l}^{\tau_{l+1}} d\tau\ \frac{2(1-\tau)}{(1+\tau)^{3-\nu_a}}T_n({\textstyle\frac{2\tau-(\tau_l+\tau_{l+1})}{\tau_{l+1}-\tau_l}})
+    \cdot\\[0.5cm]\hfill\cdot
+    \mathds 1_{\alpha_-(|\tilde k|,\tau)<\tau_{\zeta'+1}}
+    \mathds 1_{\alpha_+(|\tilde k|,\tau)>\tau_{\zeta'}}
+    \int_{\max(\tau_{l'},\alpha_-(|k|,\tau))}^{\min(\tau_{l'+1},\alpha_+(|k|,\tau))}d\sigma\ \frac{2(1-\sigma)}{(1+\sigma)^{3-\nu_b}}T_m({\textstyle\frac{2\sigma-(\tau_{l'}+\tau_{l'+1})}{\tau_{l'+1}-\tau_{l'}}})
+  \end{largearray}
+\end{equation}
+(we need the indicator functions to ensure that the bounds of the integral are correct).
+Note that $A$ is independent of $a$ and $b$, and can be computed once and for all at the beginning of execution for all values of $k_{l,j}$ (\ref{klj}).
+We then compute the integrals using Gauss-Legendre quadratures (as is proved in the next paragraph, the integrand is non singular provided $\nu_a,\nu_b\geqslant 2$).
+\bigskip
+
+\subpoint
+Note that these integrals are not singular as long as $\nu_a,\nu_b\geqslant 2$: indeed (since the only possible problems occur at $-1$, it suffices to consider the case with only one spline), $\alpha_-,\alpha_+>-1$ for $\tau>-1$, and
+\begin{equation}
+  \alpha_\pm(k,\tau)=-1+(1+\tau)\pm\frac k2(1+\tau)^2+O(1+\tau)^3
+\end{equation}
+and
+\begin{equation}
+  \int_{\alpha_-(|k|,\tau)}^{\alpha_+(|k|,\tau)}d\sigma\ \left|\frac{2(1-\sigma)}{(1+\sigma)^{3-\nu_b}}T_m(\sigma)\right|
+  \leqslant
+  4\int_{\alpha_-(|k|,\tau)}^{\alpha_+(|k|,\tau)}d\sigma\ \frac1{(1+\sigma)^{3-\nu_b}}
+\end{equation}
+which, if $\nu_b=2$, yields
+\begin{equation}
+  4\log\left(\frac{1+\alpha_+(|k|,\tau)}{1+\alpha_-(|k|,\tau)}\right)
+  =4|k|(1+\tau)+O(1+\tau)^2
+\end{equation}
+and if $\nu_b>2$, yields
+\begin{equation}
+  \frac4{\nu_b-2}\left((1+\alpha_+(|k|,\tau))^{\nu_b-2}-(1+\alpha_-(|k|,\tau))^{\nu_b-2}\right)
+  =
+  \frac4{\nu_b-2}|k|(1+\tau)^{\nu_b-1}(1+O(1+\tau))
+  .
+\end{equation}
+Therefore,
+\begin{equation}
+  \left|\frac{2(1-\tau)}{(1+\tau)^{3-\nu_a}}T_n(\tau)\int_{\alpha_-(|k|,\tau)}^{\alpha_+(|k|,\tau)}d\sigma\ \frac{2(1-\sigma)}{(1+\sigma)^{3-\nu_b}}T_m(\sigma)\right|
+  \leqslant
+  \frac8{c_{\nu_b}}
+  |k|(1-\tau)(1+\tau)^{\nu_a+\nu_b-4}(1+O(1+\tau))
+\end{equation}
+where $c_{2}:=1$ and $c_{\nu_b}:=\nu_b-2$ if $\nu_b>2$.
+The integrand is, therefore, not singular as long as $\nu_a,\nu_b\geqslant 2$.
+\bigskip
+
+\subpoint
+Evaluating convolutions at $k=0$ is not immediate, as the formula for $A(0)$ involves a bit of a computation.
+To compute $A(0)$, we expand
+\begin{equation}
+  \alpha_-(|k|,\tau)=\tau-\frac{|k|(1+\tau)^2}2+O(k^2)
+  ,\quad
+  \alpha_+(|k|,\tau)=\tau+\frac{|k|(1+\tau)^2}2+O(k^2)
+\end{equation}
+so
+\begin{equation}
+  A_{\zeta,n;\zeta',m}^{(\nu_a,\nu_b)}(0)=\mathds 1_{\zeta=\zeta'}
+  \frac1{\pi^2}
+  \int_{\tau_\zeta}^{\tau_{\zeta+1}} d\tau\ \frac{(1-\tau)^2}{(1+\tau)^{4-\nu_a-\nu_b}}T_n({\textstyle\frac{2\tau-(\tau_\zeta+\tau_{\zeta+1})}{\tau_{\zeta+1}-\tau_\zeta}})
+  T_m({\textstyle\frac{2\tau-(\tau_{\zeta}+\tau_{\zeta+1})}{\tau_{\zeta+1}-\tau_{\zeta}}})
+\end{equation}
+which we then compute using a Gauss-Legendre quadrature.
+We can then evaluate convolutions at $k=0$:
+\begin{equation}
+  (a\hat\ast b)(0)\approx\left<\mathfrak a\mathfrak b\right>
+  :=
+  \frac1{4\pi^2|k_{l,j}|}
+  \sum_{n,m=0}^P
+  \sum_{l,l'=0}^{J-1}
+  \mathfrak F_{l,n}^{(\nu_a)}(\mathfrak a)\mathfrak F_{l',m}^{(\nu_b)}(\mathfrak b)
+  A_{l,n;l',m}^{(\nu_a,\nu_b)}(0)
+  \label{avg}
+\end{equation}
+\bigskip
+
+
+\subpoint
+Let us now compute some choices of $\mathfrak a,\mathfrak b$ more explicitly.
+\medskip
+
+\subsubpoint
+Let us start with $\mathds 1_{l,n}\hat\ast a$ where $\mathds 1_{l,n}$ is the vector which has 0's everywhere except at position $(l,n)$.
+We have
+\begin{equation}
+  \mathfrak F_{l',m}^{(\nu_1)}(\mathds 1_{l,n})=\delta_{l',l}\frac{2-\delta_{m,0}}2w_n\frac{\cos({\textstyle\frac{m\pi(1+x_n)}2})}{(1-\frac{\tau_{l+1}-\tau_l}2\sin(\frac{\pi x_j}2)+\frac{\tau_{l+1}+\tau_l}2)^{\nu_1}}
+\end{equation}
+so
+\begin{equation}
+  (\mathds 1_{n,l''}\odot\mathfrak a)_{m,l}=
+  \sum_{l,p}\sum_{l'}\frac{2-\delta_{l,0}}2w_n\frac{\cos({\textstyle\frac{l\pi(1+x_n)}2})}{(1-\frac{\tau_{l''+1}-\tau_{l''}}2\sin(\frac{\pi x_n}2)+\frac{\tau_{l''+1}+\tau_{l''}}2)^{\nu_1}}A_{l'',l;l',p}^{(\nu_1,\nu_a)}(k_{l,m})\mathfrak F_{l',p}^{(\nu_a)}(\mathfrak a)
+  .
+\end{equation}
+\bigskip
+
+\subsubpoint
+We now turn to $a\hat\ast\hat v$.
+Let
+\begin{equation}
+  \Upsilon(t,|k|):=
+  \int_{||k|-t|}^{|k|+t}ds\ \frac{s\hat v(s)}{|k|}
+  \label{Upsilon}
+\end{equation}
+We have
+\begin{equation}
+  a\hat\ast\hat v=\frac 1{4\pi^2}\int_0^\infty dt\ ta(t)\Upsilon(t,|k|)
+  .
+  \label{avUpsilon}
+\end{equation}
+Following the integration scheme\-~(\ref{anyeq_integration}),
+\begin{equation}
+  (\mathfrak a\odot\mathfrak v)_{l,i}=
+  \frac1{32\pi}\sum_{l'=0}^{J-1}(\tau_{l'+1}-\tau_{l'})\sum_{j=1}^Nw_j \cos({\textstyle\frac{\pi x_j}2})(1+\vartheta_{l'}(x_j))^2\vartheta_{l'}(x_j)\mathfrak a_{l',j}\Upsilon(\vartheta(x_j),k_{l,i})
+  .
+  \label{adotv}
+\end{equation}
+At $k=0$,
+\begin{equation}
+  \left<\mathfrak a\mathfrak v\right>=
+  \frac1{32\pi}\sum_{l'=0}^{J-1}(\tau_{l'+1}-\tau_{l'})\sum_{j=1}^Nw_j \cos({\textstyle\frac{\pi x_j}2})(1+\vartheta_{l'}(x_j))^2\vartheta_{l'}(x_j)\mathfrak a_{l',j}\Upsilon(\vartheta(x_j),0)
+  \label{<av>}
+\end{equation}
+with
+\begin{equation}
+  \Upsilon(t,0)=2tv(t)
+  .
+\end{equation}
+The quantities $\Upsilon(\vartheta(x_n),k_{l,i})$ and $\Upsilon(\vartheta(x_n),0)$ are independent of $a$ and can be computed once and for all at execution.
+The integral in\-~(\ref{Upsilon}) is computed using Gauss-Legendre quadratures, but without splitting into splines.
+To maintain a high precision, we set the order of the integration to \refname{param:anyeq_J}{{\tt J}}$\times$\refname{param:anyeq_N}{{\tt N}}.
+\bigskip
+
+\subsubpoint
+Finally,
+\begin{equation}
+  (\mathds 1_{l',n}\odot\mathfrak v)_{l,i}=
+  \frac{\tau_{l'+1}-\tau_{l'}}{32\pi}w_n\cos({\textstyle\frac{\pi x_n}2})(1+\vartheta_{l'}(x_n))^2\vartheta_{l'}(x_n)\Upsilon(\vartheta(x_n),k_{l,i})
+\end{equation}
+and
+\begin{equation}
+  \left<\mathds 1_{l',n}\mathfrak v\right>=
+  \frac{\tau_{l'+1}-\tau_{l'}}{32\pi}w_n\cos({\textstyle\frac{\pi x_n}2})(1+\vartheta_{l'}(x_n))^2\vartheta_{l'}(x_n)\Upsilon(\vartheta(x_n),0)
+  .
+\end{equation}
+\bigskip
+
+\point{\bf Evaluating $\hat l_3$.}
+The only term in\-~(\ref{anyeq}) that does not involve just convolutions (whose computation was described above) is $\hat l_3$ (\ref{l3}).
+To evaluate it, we first change variables, using a generalization of bipolar coordinates (see lemma\-~\ref{lemma:bipolar2}):
+\begin{equation}
+  \begin{largearray}
+    \hat l_3(|k|)
+    =
+    \frac1{(2\pi)^5\rho^2| k|^2}\int_0^\infty dt\int_{|| k|-t|}^{| k|+t}ds\int_0^\infty dt'\int_{|| k|-t'|}^{| k|+t'}ds'\int_0^{2\pi}d\theta\ sts't'U(s)U(t)U(s')U(t')
+    \cdot\\\hfill\cdot
+    S(\mathfrak R(s,t,s',t',\theta,| k|))
+  \end{largearray}
+\end{equation}
+with
+\begin{equation}
+  \begin{largearray}
+    \mathfrak R(s,t,s',t',\theta,| k|)
+    :=
+    \frac1{\sqrt 2| k|}
+    \Big(
+      | k|^2(s^2+t^2+(s')^2+(t')^2)-| k|^4-(s^2-t^2)((s')^2-(t')^2)
+      \\[0.5cm]\hfill
+      -
+      \sqrt{(4| k|^2s^2-(| k|^2+s^2-t^2)^2)(4| k|^2(s')^2-(| k|^2+(s')^2-(t')^2)^2)}\cos\theta
+    \Big)^{\frac12}
+    .
+  \end{largearray}
+\end{equation}
+We change variables as in\-~(\ref{convvar}):
+\begin{equation}
+  \begin{largearray}
+    \hat l_3(|k|)
+    =
+    \frac1{(2\pi)^5\rho^2| k|^2}
+    \int_{-1}^1 d\tau\ \frac{2(1-\tau)}{(1+\tau)^3}U({\textstyle\frac{1-\tau}{1+\tau}})
+    \int_{\alpha_-(| k|,\tau)}^{\alpha_+(| k|,\tau)} d\sigma\ \frac{2(1-\sigma)}{(1+\sigma)^3}U({\textstyle\frac{1-\sigma}{1+\sigma}})
+    \int_{-1}^1 d\tau'\ \frac{2(1-\tau')}{(1+\tau')^3}
+    \cdot\\[0.5cm]\hfill\cdot
+    U({\textstyle\frac{1-\tau'}{1+\tau'}})
+    \int_{\alpha_-(| k|,\tau')}^{\alpha_+(| k|,\tau')} d\sigma'\ \frac{2(1-\sigma')}{(1+\sigma')^3}U({\textstyle\frac{1-\sigma'}{1+\sigma'}})
+    \int_0^{2\pi}d\theta\  S(\mathfrak R({\textstyle \frac{1-\sigma}{1+\sigma},\frac{1-\tau}{1+\tau},\frac{1-\sigma'}{1+\sigma'},\frac{1-\tau'}{1+\tau'},\theta,| k|}))
+    .
+  \end{largearray}
+\end{equation}
+We expand $U$ and $ S$ into Chebyshev polynomials as in\-~(\ref{a_chebyshev_spline}), and split the integrals into splines:
+\begin{equation}
+  \hat l_3(|k_{l,j}|)
+  \approx
+  (\mathbb U^{\otimes4}\otimes\mathbb S)_{l,j}:=
+  \sum_{\lambda_1,\cdots,\lambda_5=0}^{J-1}
+  \sum_{n_1,\cdots,n_5=1}^\infty
+  \bar A_{\lambda_1,n_1;\cdots;\lambda_5,n_5}^{(\nu)}(|k_{l,j}|)
+  \prod_{i=1}^4
+  \left(\mathfrak F_{\lambda_i,n_i}^{(\nu)}(\mathbb U)\right)
+  \mathfrak F_{\lambda_5,n_5}^{(\nu)}(\mathbb S)
+  \label{otimes}
+\end{equation}
+with $\mathbb U_{l,j}:=U(k_{l,j})$, $\mathbb S_{l,j}:=S(k_{l,j})$ and
+\begin{equation}
+  \begin{largearray}
+    \bar A_{\lambda_1,n_1;\cdots;\lambda_5,n_5}^{(\nu)}(| k|)
+    =
+    \frac1{(2\pi)^5\rho^2| k|^2}
+    \int_{\tau_{\lambda_1}}^{\tau_{\lambda_1+1}} d\tau\ \frac{2(1-\tau)}{(1+\tau)^{3-\nu}}T_{n_1}({\textstyle\frac{2\tau-(\tau_{\lambda_1}+\tau_{\lambda_1+1})}{\tau_{\lambda_1+1}-\tau_{\lambda_1}}})
+    \cdot\\[0.5cm]\hskip1cm\cdot
+    \mathds 1_{\alpha_-(| k|,\tau)<\tau_{\lambda_2+1}}\mathds 1_{\alpha_+(| k|,\tau)>\tau_{\lambda_2}}
+    \int_{\max(\tau_{\lambda_2},\alpha_-(| k|,\tau))}^{\min(\tau_{\lambda_2+1},\alpha_+(| k|,\tau))} d\sigma\ \frac{2(1-\sigma)}{(1+\sigma)^{3-\nu}}T_{n_2}({\textstyle\frac{2\sigma-(\tau_{\lambda_2}+\tau_{\lambda_2+1})}{\tau_{\lambda_2+1}-\tau_{\lambda_2}}})
+    \cdot\\[0.5cm]\hskip1cm\cdot
+    \int_{\tau_{\lambda_3}}^{\tau_{\lambda_3+1}} d\tau'\ \frac{2(1-\tau')}{(1+\tau')^{3-\nu}}T_{n_3}({\textstyle\frac{2\tau'-(\tau_{\lambda_3}+\tau_{\lambda_3+1})}{\tau_{\lambda_3+1}-\tau_{\lambda_3}}})
+    \cdot\\[0.5cm]\hskip1cm\cdot
+    \mathds 1_{\alpha_-(| k|,\tau')<\tau_{\lambda_4+1}}\mathds 1_{\alpha_+(| k|,\tau')>\tau_{\lambda_4}}
+    \int_{\max(\tau_{\lambda_4},\alpha_-(| k|,\tau'))}^{\min(\tau_{\lambda_4+1},\alpha_+(| k|,\tau'))} d\sigma'\ \frac{2(1-\sigma')}{(1+\sigma')^{3-\nu}}T_{n_4}({\textstyle\frac{2\sigma'-(\tau_{\lambda_4}+\tau_{\lambda_4+1})}{\tau_{\lambda_4+1}-\tau_{\lambda_4}}})
+    \cdot\\[0.5cm]\hskip1cm\cdot
+    \int_0^{2\pi}d\theta\ 
+    \frac{2^\nu}{(2+\mathfrak R)^\nu}
+    T_{n_5}({\textstyle\frac{2\frac{1-\mathfrak R}{1+\mathfrak R}-(\tau_{\lambda_5}+\tau_{\lambda_5+1})}{\tau_{\lambda_5+1}-\tau_{\lambda_5}}})
+    \mathds 1_{\tau_{\lambda_5}<\frac{1-\mathfrak R}{1+\mathfrak R}<\tau_{\lambda_5+1}}
+  \end{largearray}
+\end{equation}
+in which $\mathfrak R$ is a shorthand for $\mathfrak R({\textstyle \frac{1-\sigma}{1+\sigma},\frac{1-\tau}{1+\tau},\frac{1-\sigma'}{1+\sigma'},\frac{1-\tau'}{1+\tau'},\theta,| k|})$.
+Note that $\bar A$ is independent of $U$, and can be computed once and for all at the beginning of execution.
+Since the tensor $\bar A$ is quite large (it contains $(NJ)^5$ entries), and its computation can be rather long it can be inconvenient to compute $\bar A$ at every execution.
+Instead, one can use the \refname{command:anyeq_save_Abar}{{\tt save\_Abar}} method to compute $\bar A$ and save it to a file, which can then be recalled via the {\tt -s} option on the command line.
+\bigskip
+
+\point{\bf Main algorithm to compute $U$.}
+We are now ready to detail how $U\equiv\rho\hat u$ is computed.
+All of the observables we are interested in are approximated from the values $\mathbb U_{l,j}:=U(k_{l,j})$ (\ref{klj}).
+\bigskip
+
+\subpoint
+The equation for $(\mathbb U_{l,j})_{l\in\{0,\cdots,J-1\},j\in\{1,\cdots,N\}}$ is obtained by approximating\-~(\ref{anyeq_hatu}) according to the prescriptions detailed above:
+\begin{equation}
+  \mathbb U_{l,j}:=
+  \frac{1+\mathbb Y_{l,j}}{2(\mathbb X_{l,j}+1)}\ \Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
+  \label{bbU}
+\end{equation}
+\begin{equation}
+  \mathbb X_{l,j}:=\frac1{\mathbb L_{l,j}}\left(\mathbb K_{l,j}-\mathbb L_{l,j}+\frac{k_{l,j}^2}{2\rho}+\mathbb I_{l,j}\right)
+  ,\quad
+  \mathbb Y_{l,j}:=\frac1{\mathbb L_{l,j}}\left(\mathbb S_{l,j}-\mathbb L_{l,j}+\frac12\mathbb J_{l,j}+\mathbb G_{l,j}\right)
+\end{equation}
+\begin{equation}
+  \mathbb S_{l,j}:=
+  \mathfrak v_{l,j}-\frac1{\rho}(\mathbb U\odot\mathfrak v)_{l,j}
+  ,\quad
+  \mathfrak v_{l,j}:=\hat v(k_{l,j})
+  ,\quad
+  \mathbb E:=\hat v(0)-\frac1{\rho}\left<\mathbb U\mathfrak v\right>
+\end{equation}
+\begin{equation}
+  \mathbb I_{l,j}:=\frac1{\rho}\gamma_{L,2}(\beta_{L,2}(\mathbb U\odot(\mathbb U\mathbb S))_{l,j}+(1-\beta_{L,2})\mathbb E(\mathbb U\odot\mathbb U)_{l,j})
+\end{equation}
+\begin{equation}
+  \mathbb K_{l,j}:=\gamma_K(\beta_K\mathbb  S_{l,j}+(1-\beta_K)\mathbb E)
+  ,\quad
+  \mathbb L_{l,j}:=\gamma_{L,1}(\beta_{L,1}\mathbb  S_{l,j}+(1-\beta_{L,1})\mathbb E)
+\end{equation}
+\begin{equation}
+  \mathbb J_{l,j}:=\gamma_{L,3}(1-\beta_{L,3})\frac{\mathbb E}{\rho^2}(\mathbb U\odot\mathbb U)^2_{l,j}
+  +\gamma_{L,3}\beta_{L,3}(\mathbb U^{\otimes4}\otimes\mathbb S)_{l,j}
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \mathbb G_{l,j}:=
+    \frac2{\rho}\gamma_K\alpha_K(\mathbb U\odot((\beta_K\mathbb S+(1-\beta_K)\mathbb E)\mathbb U))_{l,j}
+    \\[0.5cm]\hfill
+    -
+    \frac1{\rho}\alpha_{L,1}(\mathbb U\odot((\beta_{L,1}\mathbb S+(1-\beta_{L,1})\mathbb E)\mathbb U^2))_{l,j}
+    +
+    \frac2{\rho}\alpha_{L,2}(\mathbb U\odot(\mathbb I\mathbb U))_{l,j}
+    -
+    \frac1{2\rho}\alpha_{L,3}(\mathbb U\odot \mathbb J)_{l,j}
+  \end{largearray}
+\end{equation}
+(see\-~(\ref{odot}), (\ref{otimes}) and\-~(\ref{avg}) for the definitions of $\odot$, $\otimes$ and $\left<\cdot\right>$).
+\bigskip
+
+\subpoint
+We rewrite\-~(\ref{bbU}) as a root finding problem:
+\begin{equation}
+  \Xi_{l,j}(\mathbb U):=
+  \mathbb U_{l,j}-
+  \frac{1+\mathbb Y_{l,j}}{2(\mathbb X_{l,j}+1)}\ \Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
+  =0
+  \label{root_anyeq}
+\end{equation}
+which we solve using the Newton algorithm, that is, we define a sequence of $\mathbb U$'s:
+\begin{equation}
+  \mathbb U^{(n+1)}=\mathbb U^{(n)}-(D\Xi(\mathbb U^{(n)}))^{-1}\Xi(\mathbb U^{(n)})
+\end{equation}
+where $D\Xi$ is the Jacobian of $\Xi$:
+\begin{equation}
+  (D\Xi)_{l,j;l',i}:=\frac{\partial\Xi_{l,j}}{\partial\mathbb U_{l',i}}
+  .
+\end{equation}
+\bigskip
+
+\subpoint
+We initialize the algorithm with the solution of the \refname{eq:easyeq_medeq}{Medium equation}, which is computed using the \refname{sec:easyeq}{{\tt easyeq}} method.
+However, \refname{sec:easyeq}{{\tt easyeq}} only computes $\hat u$ at the momenta given by the Gauss-Legendre quadratures\-~(\ref{easyeq_ki}).
+To obtain a value for $\hat u$ at $k_{l,j}$, we use a linear interpolation (code is provided for a polynomial interpolation, which has not performed as well).
+The parameters \refname{param:anyeq_tolerance}{{\tt tolerance}}, \refname{param:anyeq_maxiter}{{\tt maxiter}}, \refname{param:anyeq_minlrho_init}{{\tt minlrho\_init}}, \refname{param:anyeq_nlrho_init}{{\tt nlrho\_init}} are passed on to \refname{sec:easyeq}{easyeq} as is, and the order of the Gauss-Legendre quadratures is set to \refname{param:anyeq_J}{{\tt J}}$\times$\refname{param:anyeq_N}{{\tt N}}.
+\bigskip
+
+\subpoint
+We are left with computing the Jacobian of $\Xi$:
+\begin{equation}
+  \begin{largearray}
+    \partial_{l',i}\Xi_{l,j}\equiv
+    \frac{\partial\Xi_{l,j}}{\partial\mathbb U_{l',i}}=
+    \delta_{l,l'}\delta_{i,j}
+    +\frac1{2(\mathbb X_{l,j}+1)}\left(
+      \frac{(1+\mathbb Y_{l,j})\partial_{l',i}\mathbb X_{l,j}}{\mathbb X_{l,j}+1}
+      -
+      \partial_{l',i}\mathbb Y_{l,j}
+    \right)
+    \Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
+    \\[0.5cm]\hfill
+    +\frac{(1+\mathbb Y_{l,j})}{2(\mathbb X_{l,j}+1)^3}\left(
+      \frac{2(1+\mathbb Y_{l,j})}{\mathbb X_{l,j}+1}\partial_{l',i}\mathbb X_{l,j}
+      -\partial_{l',i}\mathbb Y_{l,j}
+    \right)\ \partial\Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
+  \end{largearray}
+\end{equation}
+with
+\begin{equation}
+  \partial_{l',i}\mathbb X_{l,j}=
+  \frac1{\mathbb L_{l,j}}\left(\partial_{l',i}\mathbb K_{l,j}-\partial_{l',i}\mathbb L_{l,j}+\partial_{l',i}\mathbb I_{l,j}\right)
+  -
+  \frac{\partial_{l',i}\mathbb L_{l,j}}{\mathbb L_{l,j}}\mathbb X_{l,j}
+\end{equation}
+\begin{equation}
+  \partial_{l',i}\mathbb Y_{l,j}=
+  \frac1{\mathbb L_{l,j}}\left(\partial_{l',i}\mathbb S_{l,j}-\partial_{l',i}\mathbb L_{l,j}+\frac12\partial_{l',i}\mathbb J_{l,j}+\partial_{l',i}\mathbb G_{l,j}\right)
+  -
+  \frac{\partial_{l',i}\mathbb L_{l,j}}{\mathbb L_{l,j}}\mathbb Y_{l,j}
+\end{equation}
+\begin{equation}
+  \partial_{l',i}\mathbb S_{l,j}=
+  -\frac1{\rho}(\mathds 1_{l',i}\odot\mathfrak v)_{l,j}
+  ,\quad
+  \partial_{l',i}\mathbb E=-\frac1{\rho}\left<\mathds 1_{l',i}\mathfrak v\right>
+\end{equation}
+\begin{equation}
+  \partial_{l',i}\mathbb K_{l,j}:=\gamma_K(\beta_K\partial_{l',i}\mathbb  S_{l,j}+(1-\beta_K)\partial_{l',i}\mathbb E)
+  ,\quad
+  \partial_{l',i}\mathbb L_{l,j}:=\gamma_{L,1}(\beta_{L,1}\partial_{l',i}\mathbb  S_{l,j}+(1-\beta_{L,1})\partial_{l',i}\mathbb E)
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \partial_{l',i}\mathbb I_{l,j}=
+    \frac1{\rho}\gamma_{L,2}\left(\beta_{L,2}(
+      \mathds 1_{l',i}\odot(\mathbb U\mathbb S)
+      +
+      \mathbb U\odot(\mathds 1_{l',i}\mathbb S+\mathbb U\partial_{l',i}\mathbb S)
+    )_{l,j}
+    \right.\\[0.5cm]\hfill\left.
+    +
+    (1-\beta_{L,2})(
+      \partial_{l',i}\mathbb E(\mathbb U\odot\mathbb U)_{l,j}
+      +
+      2\mathbb E(\mathbb U\odot\mathds 1_{l',i})_{l,j}
+    )\right)
+  \end{largearray}
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \partial_{l',i}\mathbb J_{l,j}=
+    \frac1{\rho^2}\gamma_{L,3}(1-\beta_{L,3})\left(
+      \partial_{l',i}\mathbb E(\mathbb U\odot\mathbb U)^2_{l,j}
+      +
+      4\mathbb E(\mathbb U\odot\mathbb U)_{l,j}(\mathbb U\odot\mathds 1_{l',i})_{l,j}
+    \right)
+    \\[0.5cm]\hfill
+    +\gamma_{L,3}\beta_{L,3}(
+      4\mathds 1_{l',i}\otimes\mathbb U^{\otimes3}\otimes\mathbb S
+      +
+      \mathbb U^{\otimes4}\otimes\partial_{l',i}\mathbb S
+    )_{l,j}
+  \end{largearray}
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \partial_{l',i}\mathbb G_{l,j}=
+    \frac2{\rho}\gamma_K\alpha_K\beta_K\left(
+      \mathds 1_{l',i}\odot(\mathbb S\mathbb U)
+      +
+      \mathbb U\odot(\mathbb S\mathds 1_{l',i}
+      +\partial_{l',i}\mathbb S\mathbb U)
+    \right)_{l,j}
+    \\[0.5cm]\hfill
+    +\frac2{\rho}\gamma_K\alpha_K(1-\beta_K)\left(
+      \partial_{l',i}\mathbb E(\mathbb U\odot\mathbb U)_{l,j}
+      +2\mathbb E(\mathds 1_{l',i}\odot\mathbb U)_{l,j}
+    \right)
+    \\[0.5cm]\indent
+    -
+    \frac1{\rho}\alpha_{L,1}\beta_{L,1}\left(
+      \mathds 1_{l',i}\odot(\mathbb S\mathbb U^2)
+      +
+      \mathbb U\odot(
+	\partial_{l',i}\mathbb S\mathbb U^2
+	+
+	2\mathbb S\mathbb U\mathds 1_{l',i}
+      )
+    \right)_{l,j}
+    \\[0.5cm]\hfill
+    -
+    \frac1{\rho}\alpha_{L,1}(1-\beta_{L,1})\left(
+      \mathbb E\mathds 1_{l',i}\odot\mathbb U^2
+      +
+      \mathbb U\odot(
+	\partial_{l',i}\mathbb E\mathbb U^2
+	+
+	2\mathbb E\mathbb U\mathds 1_{l',i}
+      )
+    \right)_{l,j}
+    \\[0.5cm]\hfill
+    +
+    \frac2{\rho}\alpha_{L,2}(
+      (\mathds 1_{l',i}\odot(\mathbb I\mathbb U))_{l,j}
+      +
+      \mathbb U\odot(
+	\partial_{l',i}\mathbb I\mathbb U
+	+
+	\mathbb I\mathds 1_{l',i})
+    )_{l,j}
+    -
+    \frac1{2\rho}\alpha_{L,3}(
+      \mathds 1_{l',i}\odot \mathbb J
+      +
+      \mathbb U\odot\partial_{l',i}\mathbb J
+    )_{l,j}
+    .
+  \end{largearray}
+\end{equation}
+\bigskip
+
+\subpoint
+We iterate the Newton algorithm until the Newton relative error $\epsilon$ becomes smaller than the \refname{param:easyeq_tolerance}{{\tt tolerance}} parameter.
+The Newton error is defined as
+\begin{equation}
+  \epsilon=\frac{\|\mathbb U^{(n+1)}-\mathbb U^{(n)}\|_2}{\|\mathbb U^{(n)}\|_2}
+\end{equation}
+where $\|\cdot\|_2$ is the $l_2$ norm.
+The energy thus obtained is
+\begin{equation}
+  e=\frac\rho2\mathds E
+\end{equation}
+the Fourier transform $\hat u$ of the solution is
+\begin{equation}
+  \hat u(k_{l,j})\approx\frac{\mathbb U_{l,j}}\rho
+\end{equation}
+where $k_{l,j}$ was defined in\-~(\ref{klj}), and the solution $u$ in real space is obtained by inverting the Fourier transform, following the prescription of\-~(\ref{anyeq_integration}):
+\begin{equation}
+  u(|x|)=\int \frac{dk}{(2\pi)^3}\ e^{-ikx}\hat u(|k|)
+  \approx
+  \sum_{l=0}^{J-1}
+  \frac{\tau_{l+1}-\tau_l}{16\pi}
+  \sum_{j=1}^Nw_j
+  \cos({\textstyle\frac{\pi x_j}2})(1+\vartheta_l(x_j))^2\vartheta_l(x_j)\mathbb U_{l,j}\sin(\vartheta_l(x_j)|x|)
+  .
+\end{equation}
+\bigskip
+
+\point{\bf Condensate fraction.}
+To compute the condensate fraction, we solve the modified {\tt anyeq} (see\-~\cite{CJL20b}):
+\begin{equation}
+  (-\Delta+2\mu)u_\mu
+  =
+  (1-u_\mu)v-2\rho K+\rho^2 L
+  .
+\end{equation}
+where $K$ and $L$ are defined as in\-~(\ref{anyeq_K})-(\ref{anyeq_e}) in which $u$ is replaced with $u_\mu$.
+The uncondensed fraction is then
+\begin{equation}
+  \eta=\partial_\mu e|_{\mu=0}
+  =-\frac\rho2\int dx\ v(x)\partial_\mu u_\mu(x)|_{\mu=0}
+  .
+\end{equation}
+To compute the energy in the presence of the parameter $\mu$, we proceed in the same way as for $\mu=0$, the only difference being that $k^2$ should formally be replaced by $k^2+2\mu$.
+In other words, we consider $\mathbb U_{j,l}=u_\mu(|k_{j,l}|)$ and define $\Xi(\mathbb U,\mu)$ in the same way as in\-~(\ref{root_anyeq}), except that $\mathbb X_i$ should be replaced by
+\begin{equation}
+  \mathbb X_{l,j}=\frac1{\mathbb L_{l,j}}\left(\mathbb K_{l,j}-\mathbb L_{l,j}+\frac{k_{l,j}^2+2\mu}{2\rho}+\mathbb I_{l,j}\right)
+  .
+\end{equation}
+We then solve
+\begin{equation}
+  \Xi(\mathbb U_\mu,\mu)=0
+\end{equation}
+By differentiating this identity with respect to $\mu$, we find $\partial_\mu u_\mu$:
+\begin{equation}
+  \partial_\mu \mathbb U|_{\mu=0}=-(D\Xi)^{-1}\partial_\mu\Xi|_{\mu=0}
+\end{equation}
+and
+\begin{equation}
+  \partial_\mu\Xi|_{\mu=0}=
+  \frac{(1+\mathbb Y_{l,j})\partial_\mu\mathbb X_{l,j}}{2(\mathbb X_{l,j}+1)^2}
+  \Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
+  +\frac{(1+\mathbb Y_{l,j})^2}{(\mathbb X_{l,j}+1)^4}\partial_\mu\mathbb X_{l,j}
+  \partial\Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
+  ,\quad
+  \partial_\mu\mathbb X_{l,j}=\frac1{\rho\mathbb L_{l,j}}
+  .
+\end{equation}
+We then approximate
+\begin{equation}
+  \eta\approx
+  -\frac12\left<\mathfrak v\partial_\mu\mathbb U\right>
+\end{equation}
+(see\-~(\ref{avg})).
+\bigskip
+
+\point{\bf Correlation function.}
+The two-point correlation function is
+\begin{equation}
+  C_2(x):=
+  2\rho\frac{\delta e}{\delta v(x)}
+  .
+\end{equation}
+In Fourier space,
+\begin{equation}
+  C_2(x)=
+  2\rho\int dk\ e^{ikx}\frac{\delta e}{\delta\hat v(k)}
+  =2\rho\int dk\ \frac{\sin(|k||x|)}{|k||x|}\frac{\delta e}{\delta\hat v(k)}
+  .
+  \label{C2fourier}
+\end{equation}
+\bigskip
+
+\subpoint
+We can compute this quantity by considering a modified {\tt anyeq} in Fourier space, by formally replacing $\hat v$ with
+\begin{equation}
+  \hat v+\lambda g(|k|)
+  ,\quad
+  g(|k|):=\frac{\sin(|k||x|)}{|k||x|}
+  .
+\end{equation}
+Indeed, if $e_\lambda$ denotes the energy of this modified equation,
+\begin{equation}
+  \partial_\lambda e_\lambda|_{\lambda=0}
+  =\int dk\ \frac{\delta e}{\hat v(k)}\partial_\lambda({\textstyle \hat v(k)+\lambda g(|k|)})
+  =\int dk\ g(|k|)\frac{\delta e}{\hat v(k)}
+\end{equation}
+so, denoting the solution of the modified equation by $u_\lambda$,
+\begin{equation}
+  C_2(x)=2\rho\partial_\lambda e_\lambda|_{\lambda=0}
+  =-\rho^2\int dx\ (g(x)u_0(x)+v(x)\partial_\lambda u_\lambda(x)|_{\lambda=0})
+  .
+\end{equation}
+We compute $\partial_\lambda u_\lambda|_{\lambda=0}$ in the same way as the uncondensed fraction: we define $\Xi(\mathbb U,\lambda)$ by formally adding $\lambda g(|k|)$ to $\hat v$, solve $\Xi(\mathbb U,\lambda)=0$, and differentiate:
+\begin{equation}
+  \partial_\lambda\mathbb U|_{\lambda=0}=-(D\Xi)^{-1}\partial_\lambda\Xi|_{\lambda=0}
+  .
+\end{equation}
+\bigskip
+
+\subpoint
+We compute $\partial_\lambda\Xi|_{\lambda=0}$:
+\begin{equation}
+  \begin{largearray}
+    \partial_\lambda\Xi_{l,j}=
+    \delta_{l,l'}\delta_{i,j}
+    +\frac1{2(\mathbb X_{l,j}+1)}\left(
+      \frac{(1+\mathbb Y_{l,j})\partial_\lambda\mathbb X_{l,j}}{\mathbb X_{l,j}+1}
+      -
+      \partial_\lambda\mathbb Y_{l,j}
+    \right)
+    \Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
+    \\[0.5cm]\hfill
+    +\frac{(1+\mathbb Y_{l,j})}{2(\mathbb X_{l,j}+1)^3}\left(
+      \frac{2(1+\mathbb Y_{l,j})}{\mathbb X_{l,j}+1}\partial_\lambda\mathbb X_{l,j}
+      -\partial_\lambda\mathbb Y_{l,j}
+    \right)\ \partial\Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
+  \end{largearray}
+\end{equation}
+with
+\begin{equation}
+  \partial_\lambda\mathbb S_{l,j}=
+  \mathfrak g_{l,j}
+  -\frac1{\rho}(\mathbb U\odot\mathfrak g)_{l,j}
+  ,\quad
+  \partial_\lambda\mathbb E=g(0)-\frac1{\rho}\left<\mathbb U\mathfrak g\right>
+\end{equation}
+\begin{equation}
+  \partial_\lambda\mathbb X_{l,j}=
+  \frac1{\mathbb L_{l,j}}\left(\partial_\lambda\mathbb K_{l,j}-\partial_\lambda\mathbb L_{l,j}+\partial_\lambda\mathbb I_{l,j}\right)
+  -
+  \frac{\partial_\lambda\mathbb L_{l,j}}{\mathbb L_{l,j}}\mathbb X_{l,j}
+\end{equation}
+\begin{equation}
+  \partial_\lambda\mathbb Y_{l,j}=
+  \frac1{\mathbb L_{l,j}}\left(\partial_\lambda\mathbb S_{l,j}-\partial_\lambda\mathbb L_{l,j}+\frac12\partial_\lambda\mathbb J_{l,j}+\partial_\lambda\mathbb G_{l,j}\right)
+  -
+  \frac{\partial_\lambda\mathbb L_{l,j}}{\mathbb L_{l,j}}\mathbb Y_{l,j}
+\end{equation}
+\begin{equation}
+  \partial_\lambda\mathbb K_{l,j}:=\gamma_K(\beta_K\partial_\lambda\mathbb  S_{l,j}+(1-\beta_K)\partial_\lambda\mathbb E)
+  ,\quad
+  \partial_\lambda\mathbb L_{l,j}:=\gamma_{L,1}(\beta_{L,1}\partial_\lambda\mathbb  S_{l,j}+(1-\beta_{L,1})\partial_\lambda\mathbb E)
+\end{equation}
+\begin{equation}
+  \partial_\lambda\mathbb I_{l,j}=
+  \frac1{\rho}\gamma_{L,2}\left(
+    \beta_{L,2}(\mathbb U\odot(\mathbb U\partial_\lambda\mathbb S)_{l,j}
+    +
+    (1-\beta_{L,2})\partial_\lambda\mathbb E(\mathbb U\odot\mathbb U)_{l,j}
+  \right)
+\end{equation}
+\begin{equation}
+  \partial_\lambda\mathbb J_{l,j}=
+  \frac1{\rho^2}\gamma_{L,3}(1-\beta_{L,3})\partial_\lambda\mathbb E(\mathbb U\odot\mathbb U)^2_{l,j}
+  +\gamma_{L,3}\beta_{L,3}(\mathbb U^{\otimes4}\otimes\partial_\lambda\mathbb S)_{l,j}
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \partial_\lambda\mathbb G_{l,j}=
+    \frac2{\rho}\gamma_K\alpha_K\beta_K\left(\partial_\lambda\mathbb S\mathbb U\right)_{l,j}
+    +\frac2{\rho}\gamma_K\alpha_K(1-\beta_K)\partial_\lambda\mathbb E(\mathbb U\odot\mathbb U)_{l,j}
+    -
+    \frac1{\rho}\alpha_{L,1}\beta_{L,1}\left(\mathbb U\odot\partial_\lambda\mathbb S\mathbb U^2\right)_{l,j}
+    \\[0.5cm]\hfill
+    -
+    \frac1{\rho}\alpha_{L,1}(1-\beta_{L,1})\partial_\lambda\mathbb E\left(\mathbb U\odot(\mathbb U^2)\right)_{l,j}
+    +
+    \frac2{\rho}\alpha_{L,2}\mathbb U\odot(\partial_\lambda\mathbb I\mathbb U)_{l,j}
+    -
+    \frac1{2\rho}\alpha_{L,3}(\mathbb U\odot\partial_\lambda\mathbb J)_{l,j}
+    .
+  \end{largearray}
+\end{equation}
+To evaluate $(\mathbb U\odot\mathfrak g)$ and $\left<\mathbb U\mathfrak g\right>$, we proceed as in\-~(\ref{adotv}) and\-~(\ref{<av>}).
+To do so, we replace $\hat v$ with $g$ in the computation of $\Upsilon$.
+\bigskip
+
+\subpoint
+In order to invert the Fourier transform in\-~(\ref{C2fourier}) numerically, we will use a Hann window (see appendix\-~\ref{appendix:hann})
+\begin{equation}
+  H_L(k):=\mathds 1_{|k|<\frac L2}\cos^2({\textstyle\frac{\pi|k|}{L}})
+  .
+  \label{hann}
+\end{equation}
+The parameter $L$ is set using \refname{param:anyeq_window_L}{{\tt window\_L}}.
+The computation is changed only in that $g$ is changed to $H_L(k)\frac{\sin(|k||x|)}{|k||x|}$.
+\bigskip
+
+\point {\bf Momentum distribution.}
+To compute the momentum distribution (see\-~\cite{CHe21}), we add a parameter $\lambda$ to {\tt anyeq}:
+\begin{equation}
+  -\Delta u_\lambda(|x|)
+  =
+  (1-u_\lambda(|x|))v(|x|)-2\rho K(|x|)+\rho^2 L(|x|)
+  -2\lambda \hat u_0(q)\cos(q\cdot x)
+\end{equation}
+($\hat u_0\equiv\hat u_\lambda|_{\lambda=0}$).
+The momentum distribution is then
+\begin{equation}
+  \mathcal M(q)=\partial_\lambda e|_{\lambda=0}
+  =-\frac\rho2\int dx\ v(x)\partial_\lambda u_\lambda(x)|_{\lambda=0}
+  .
+\end{equation}
+Note that the Fourier transform of $2\lambda\hat u_0(q)\cos(q\cdot x)$ is
+\begin{equation}
+  -(2\pi)^3\lambda\hat u_0(q)(\delta(q+k)+\delta(q-k))
+  .
+  \label{fouriermomentum}
+\end{equation}
+We compute $\partial_\lambda u_\lambda|_{\lambda=0}$ in the same way as the uncondensed fraction.
+\bigskip
+
+\subpoint
+In order to do so we will discretize momentum space, see\-~(\ref{klj}), and so it is necessary to construct a discrete analog of the delta-functions in\-~(\ref{fouriermomentum}).
+The starting point we take is
+\begin{equation}
+  \int dk f(k)\delta(k-q)
+  =f(q)
+\end{equation}
+so, when approximating the integral according to\-~(\ref{anyeq_integration}), we find
+\begin{equation}
+  \frac{\pi^2}2\sum_{l=0}^{J-1}
+  (\tau_{l+1}-\tau_l)
+  \sum_{j=1}^Nw_j
+  \cos({\textstyle\frac{\pi x_j}2})(1+\vartheta_l(x_j))^2\vartheta_l^2(x_j)f(\vartheta_l(x_j))\tilde\delta(\vartheta_l(x_j)-k_{l',i})
+  \approx
+  f(k_{l',i})
+\end{equation}
+where $\tilde\delta$ is the approximation of the delta-function.
+Since
+\begin{equation}
+  \vartheta_l(x_j)\equiv k_{l,j}
+\end{equation}
+(see\-~(\ref{chebyshevvars})), we find that the definition of $\tilde\delta$ must be
+\begin{equation}
+  \tilde\delta_{l,j;l',i}
+  :=
+  \delta_{l,l'}\delta_{j,i}
+  \frac2{\pi^2}
+  \left(
+    (\tau_{l+1}-\tau_l)
+    w_j
+    \cos({\textstyle\frac{\pi x_j}2})
+    (1+k_{l,j})^2k_{l,j}^2
+  \right)^{-1}
+  .
+\end{equation}
+Note that, due to the invariance under rotation, the approximation for $\delta(q+k)$ is equal to that for $\delta(q-k)$.
+\bigskip
+
+\subpoint
+To compute the momentum distribution at $q=k_{l',i}$, we define $\Xi(\mathbb U,\lambda)$ by formally adding $-2(2\pi)^3\lambda \hat u_0(k_{l',i})\tilde\delta_{l,j;l',i}$ to $\mathbb G_{l,j}$, which corresponds to replacing $\mathbb Y$ with
+\begin{equation}
+  \mathbb Y_{l,j}=
+  \frac1{\mathbb L_{l,j}}\left(\mathbb S_{l,j}-\mathbb L_{l,j}+\frac12\mathbb J_{l,j}+\mathbb G_{l,j}-2(2\pi)^3\lambda\frac{\mathbb U_{l',i}}\rho\tilde\delta_{l,j;l',i}\right)
+\end{equation}
+Then we solve $\Xi(\mathbb U,\lambda)=0$, and differentiate:
+\begin{equation}
+  \partial_\lambda\mathbb U|_{\lambda=0}=-(D\Xi)^{-1}\partial_\lambda\Xi|_{\lambda=0}
+  .
+\end{equation}
+Finally,
+\begin{equation}
+  \partial_\lambda\Xi_{l,j}|_{\lambda=0}
+  =
+  -\partial_\lambda\mathbb Y_{l,j}|_{\lambda=0}\left(
+    \frac1{2(\mathbb X_{l,j}+1)}
+    \Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
+    +\frac{(1+\mathbb Y_{l,j})}{2(\mathbb X_{l,j}+1)^3}
+    \partial\Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
+  \right)
+\end{equation}
+with
+\begin{equation}
+  \partial_\lambda\mathbb Y_{l,j}|_{\lambda=0}=
+  \frac{2(2\pi)^3\mathbb U_{l',i}\tilde\delta_{l,j;l',i}}{\rho\mathbb L_{l,j}}
+  .
+\end{equation}
+
+\subsection{{\tt simpleq-Kv}}\label{sec:simpleq-Kv}
+\indent
+The method is used to compute observables for the simple equation
+\begin{equation}
+  -\Delta u
+  =
+  v(1-u)-4eu+2e\rho u\ast u
+  ,\quad
+  e=\frac\rho2\int dx\ (1-u(|x|))v(|x|)
+  .
+  \label{simpleq}
+\end{equation}
+One can show\-~\cite[Theorem\-~1.6]{CJL20b} that the condensate fraction is
+\begin{equation}
+  \eta=\frac{\rho\int v(x)\mathfrak Ku(x)\ dx}{1-\rho\int v(x)\mathfrak K(2u(x)-\rho u\ast u(x))\ dx}
+  \label{simpleq-Kv_eta}
+\end{equation}
+with
+\begin{equation}\label{bg6}
+  \mathfrak K  = (-\Delta + v +  4e(1 - \rho u\ast))^{-1}
+  .
+  \label{Kv}
+\end{equation}
+Similarly, the two-point correlation function is\-~\cite[(45)]{CHe21}
+\begin{equation}
+  C_2=
+  \rho^2(1-u)+
+  \rho^2\frac{\mathfrak Kv(1-u)-2\rho u\ast \mathfrak Kv+\rho^2u\ast u\ast \mathfrak Kv}{1-\rho\int dx\ v(x)\mathfrak K(2u(x)-\rho u\ast u(x))}
+  .
+  \label{simpleq-Kv_C2}
+\end{equation}
+Thus, using the fact that $\mathfrak K$ is self-adjoint, we can compute these observables of the simple equation directly from the knowledge of $\mathfrak Kv$.
+\bigskip
+
+\subsubsection{Usage}
+\indent
+The computation uses the same approximation scheme as \refname{sec:anyeq}{{\tt anyeq}}, as well as using the solution of \refname{sec:anyeq}{{\tt anyeq}}.
+As such, it takes a similar list of parameters:
+\refname{param:anyeq_rho}{{\tt rho}},
+\refname{param:anyeq_tolerance}{{\tt tolerance}},
+\refname{param:anyeq_maxiter}{{\tt maxiter}},
+\refname{param:anyeq_P}{{\tt P}},
+\refname{param:anyeq_N}{{\tt N}},
+\refname{param:anyeq_J}{{\tt J}},
+\refname{param:anyeq_minlrho_init}{{\tt minlrho\_init}},
+\refname{param:anyeq_nlrho_init}{{\tt nlrho\_init}}.
+\bigskip
+
+The available {\tt commands} are the following.
+\bigskip
+
+\begin{itemize}
+  \item\makelink{command:simpleq-Kv_Kv}{}
+  {\tt Kv}: compute $\mathfrak Kv$ as a function of $|x|$.\par
+  \underline{Extra parameters}:
+  \begin{itemize}
+    \item\makelink{param:simpleq-Kv_xmin}{}
+    {\tt xmin} ({\tt Float64}, default: 0): minimum of the range of $|x|$ to be printed.
+
+    \item\makelink{param:simpleq-Kv_xmax}{}
+    {\tt xmax} ({\tt Float64}, default: 100): maximum of the range of $|x|$ to be printed.
+
+    \item\makelink{param:simpleq-Kv_nx}{}
+    {\tt nx} ({\tt Int64}, default: 100): number of points to print (linearly spaced).
+
+  \end{itemize}
+  \underline{Output} (one line for each value of $|x|$): [$|x|$] [$\mathfrak K v$]
+
+  \item\makelink{command:simpleq-Kv_condensate-fraction}{}
+  {\tt condensate-fraction}: compute the uncondensed fraction as a function of $\rho$\par
+  \underline{Extra parameters}:
+  \begin{itemize}
+    \item\makelink{param:simpleq-Kv_minlrho}{}
+    {\tt minlrho} ({\tt Float64}, default: $10^{-6}$): minimal value for $\log_{10}\rho$.
+
+    \item\makelink{param:simpleq-Kv_maxlrho}{}
+    {\tt maxlrho} ({\tt Float64}, default: $10^{2}$): maximal value for $\log_{10}\rho$.
+
+    \item\makelink{param:simpleq-Kv_nlrho}{}
+    {\tt nlrho} ({\tt Int64}, default: $100$): number of values for $\rho$ (spaced logarithmically).
+
+    \item\makelink{param:simpleq-Kv_rhos}{}
+    {\tt rhos} ({\tt Array\{Float64\}}, default: $(10^{{\tt minlrho}+\frac{{\tt maxlrho}-{\tt minlrho}}{{\tt nlrho}}n})_n$: list of values for $\rho$, specified as a `{\tt,}' separated list.
+  \end{itemize}
+  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$\eta$]
+
+  \item\makelink{command:simpleq-Kv_2pt}{}
+  {\tt 2pt}: compute the two-point correlation as a function of $|x|$.\par
+  \underline{Extra parameters}: same as \refname{command:simpleq-Kv_Kv}{{\tt Kv}}.\par
+  \underline{Output} (one line for each value of $|x|$): [$|x|$] [$C_2$]
+\end{itemize}
+
+\subsubsection{Description}
+\indent
+In Fourier space\-~(\ref{anyeq_fourier}),
+\begin{equation}
+  \hat{\mathfrak K}^{-1}f
+  :=
+  \int dx\ e^{ikx}
+  (-\Delta+4e(1-\rho u\ast)+v)f
+  =(k^2+4e(1-\rho\hat u(|k|))+\hat v\hat\ast)\hat f
+\end{equation}
+where $\hat\ast$ is defined in\-~(\ref{anyeq_hatast}).
+We follow the same approximation scheme as \refname{sec:anyeq}{{\tt anyeq}}:
+\begin{equation}
+  \hat{\mathfrak K}^{-1}f(k_{l,j})
+  \approx(k_{l,j}^2+4e(1-\mathbb U_{l,j}))\mathfrak f_{l,j}
+  +
+  (\mathfrak v\odot\mathfrak f)_{l,j}
+\end{equation}
+with $\mathfrak f_{l,j}:=f(k_{l,j})$, $\mathbb U_{l,j}:=\rho u(|k_{l,j}|)$, $\odot$ is defined in\-~(\ref{odot}), and $k_{l,j}$ is defined in\-~(\ref{klj}).
+Therefore, we approximate the operator $\hat{\mathfrak K}^{-1}$ by a matrix:
+\begin{equation}
+  \hat{\mathfrak K}^{-1}f(k_{l,j})\approx
+  \sum_{l'=0}^{J-1}\sum_{i=1}^N M_{l,j;l',i}\mathfrak f_{l',i}
+\end{equation}
+with, by\-~(\ref{odot}) and\-~(\ref{frakF}),
+\begin{equation}
+  \begin{largearray}
+    M_{l,j;l',i}:=
+    \delta_{l',l}\delta_{j,i}(k_{l,j}^2+4e(1-\mathbb U_{l,j}))
+    \\\hfill
+    +
+    \frac1{4\pi^2|k_{l,j}|}\sum_{n,m=0}^P\sum_{l''=0}^{J-1}\mathfrak F_{l'',n}^{(\nu_v)}(\mathfrak v)
+    A_{l'',n;l',m}^{(\nu_v,\nu_f)}(|k_{l,j}|)
+    \frac{\frac{2-\delta_{m,0}}2w_i\cos(\frac{m\pi(1+x_i)}2)}{(1-\frac{\tau_{l'+1}-\tau_{l'}}2\sin(\frac{\pi x_i}2)+\frac{\tau_{l'+1}+\tau_{l'}}2)^{\nu_f}}
+    .
+  \end{largearray}
+\end{equation}
+Defining $\mathds 1_{l',i}$ as the vector whose only non-vanishing component is that indexed by $l',i$ which is equal to 1, we can rewrite
+\begin{equation}
+  M_{l,j;l',i}:=
+  \delta_{l',l}\delta_{j,i}(k_{l,j}^2+4e(1-\mathbb U_{l,j}))
+  +
+  (\mathfrak v\odot\mathds 1_{l',i})_{l,j}
+  .
+\end{equation}
+Thus
+\begin{equation}
+  \mathfrak Kv\approx M^{-1}\mathfrak v
+  .
+\end{equation}
+\bigskip
+
+\indent
+To compute\-~(\ref{simpleq-Kv_eta}), we write
+\begin{equation}
+  \eta=\frac{\rho\int\frac{dk}{(2\pi)^3} \hat u(k)\hat{\mathfrak K}\hat v(k)}{1-\rho\int \frac{dk}{(2\pi)^3)}(2\hat u(k)-\rho\hat u^2(k))\hat{\mathfrak K}\hat v(k)}
+\end{equation}
+which we approximate as
+\begin{equation}
+  \eta\approx\frac{\left<\mathbb UM^{-1}\mathfrak v\right>}{1-\left<(2\mathbb U-\mathbb U^2)M^{-1}\mathfrak v\right>}
+\end{equation}
+where $\left<\cdot\right>$ is defined in\-~(\ref{avg}).
+We can thus compute $\eta$ using the solution $\mathbb U$ computed by \refname{sec:anyeq}{{\tt anyeq}}.
+\bigskip
+
+\indent
+To compute\-~(\ref{simpleq-Kv_C2}), we write
+\begin{equation}
+  C_2=
+  \rho^2(1-u(x))+
+  \rho^2\frac{\mathfrak Kv(x)(1-u(x))+\int\frac{dk}{(2\pi)^3}(-2\rho\hat u\hat{\mathfrak K}\hat v+\rho^2\hat u^2\hat{\mathfrak K}\hat v)}{1-\rho\int\frac{dk}{(2\pi)^3}\ (2\hat u(k)-\rho \hat u^2(k))\hat{\mathfrak K}\hat v(k)}
+\end{equation}
+which we approximate as
+\begin{equation}
+  C_2\approx
+  \rho^2-\rho\left<e^{-ik_{l,j}x}\mathbb U\right>+
+  \rho^2\frac{\left<M^{-1}\mathfrak v\right>(1-\rho^{-1}\left<e^{-ik_{l,j}x}\mathbb U\right>)+\left<(-2\mathbb U M^{-1}\mathfrak v+\mathbb UM^{-1}\mathfrak v)\right>}{1-\left<(2\mathbb U-\mathbb U^2)M^{-1}\mathfrak v\right>}
+  .
+\end{equation}
+We can thus compute $C_2$ using the solution $\mathbb U$ computed by \refname{sec:anyeq}{{\tt anyeq}}.
+
+\subsection{\tt simpleq-hardcore}\label{sec:simpleq-hardcore}
+\indent
+This method is used to solve the Simple equation with a hardcore potential:
+\begin{equation}
+  \left\{\begin{array}{>\displaystyle ll}
+    (-\Delta+4e) u(x)=2e\rho u\ast u(x)
+    &\mathrm{for\ }|x|>1\\
+    u(x)=1
+    &\mathrm{for\ }|x|\leqslant1
+  \end{array}\right.
+  \label{linearhc}
+\end{equation}
+with
+\begin{equation}
+  e=-\frac{4\pi\rho\partial u|_{|x|\searrow 1}}{2(1-\frac83\pi\rho+\rho^2\int_{|x|<1}dx\ u\ast u(x))}
+  .
+  \label{hardcore_energy}
+\end{equation}
+\bigskip
+
+This equation is solved in $x$-space, and as such is very different from \refname{sec:easyeq}{{\tt easyeq}}, and significantly longer to run.
+\bigskip
+
+\subsubsection{Usage}
+\indent
+Unless otherwise noted, this method takes the following parameters (specified via the {\tt [-p params]} flag, as a `{\tt;}' separated list of entries).
+\begin{itemize}
+  \item\makelink{param:simpleq-hardcore_rho}{}
+  {\tt rho} ({\tt Float64}, default: $10^{-6}$): density $\rho$.
+
+  \item\makelink{param:simpleq-hardcore_tolerance}{}
+  {\tt tolerance} ({\tt Float64}, default: $10^{-11}$): maximal size of final step in Newton iteration.
+
+  \item\makelink{param:simpleq-hardcore_maxiter}{}
+  {\tt maxiter} ({\tt Int64}, default: 21): maximal number of iterations of the Newton algorithm before giving up.
+
+  \item\makelink{param:simpleq-hardcore_P}{}
+  {\tt P} ({\tt Int64}, default: 11): order of all Chebyshev polynomial expansions (denoted by $P$ below).
+
+  \item\makelink{param:simpleq-hardcore_N}{}
+  {\tt N} ({\tt Int64}, default: 12): order of all Gauss quadratures (denoted by $N$ below).
+
+  \item\makelink{param:simpleq-hardcore_J}{}
+  {\tt J} ({\tt Int64}, default: 10): number of splines (denoted by $J$ below).
+\end{itemize}
+\bigskip
+
+The available {\tt commands} are the following.
+
+\begin{itemize}
+  \item\makelink{command:simpleq-hardcore_energy_rho}{}
+  {\tt energy\_rho}: compute the energy $e$ as a function of $\rho$.\par
+  \underline{Disabled parameters}: \refname{param:simpleq-hardcore_rho}{{\tt rho}}.\par
+  \underline{Extra parameters}:
+  \begin{itemize}
+    \item\makelink{param:simpleq-hardcore_minlrho}{}
+    {\tt minlrho} ({\tt Float64}, default: $10^{-6}$): minimal value for $\log_{10}\rho$.
+
+    \item\makelink{param:simpleq-hardcore_maxlrho}{}
+    {\tt maxlrho} ({\tt Float64}, default: $10^{2}$): maximal value for $\log_{10}\rho$.
+
+    \item\makelink{param:simpleq-hardcore_nlrho}{}
+    {\tt nlrho} ({\tt Int64}, default: $100$): number of values for $\rho$ (spaced logarithmically).
+
+    \item\makelink{param:simpleq-hardcore_rhos}{}
+    {\tt rhos} ({\tt Array\{Float64\}}, default: $(10^{{\tt minlrho}+\frac{{\tt maxlrho}-{\tt minlrho}}{{\tt nlrho}}n})_n$: list of values for $\rho$, specified as a `{\tt,}' separated list.
+  \end{itemize}
+  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$e$] [Newton error $\epsilon$].\par
+  \underline{Multithread support}: yes, different values of $\rho$ split up among workers.
+
+  \item\makelink{command:simpleq-hardcore_condensate_fraction_rho}{}
+  {\tt condensate\_fraction\_rho}: compute the uncondensed fraction $\eta$ as a function of $\rho$.\par
+  \underline{Disabled parameters}: same as \refname{command:simpleq-hardcore_energy_rho}{{\tt energy\_rho}}.\par
+  \underline{Extra parameters}: same as \refname{command:simpleq-hardcore_energy_rho}{{\tt energy\_rho}}.\par
+  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$\eta$] [Newton error $\epsilon$].\par
+  \underline{Multithread support}: yes, different values of $\rho$ split up among workers.
+
+  \item\makelink{command:simpleq-hardcore_ux}{}
+  {\tt ux}: compute $u$ as a function of $|x|$.\par
+  \underline{Extra parameters}:
+  \begin{itemize}
+    \item\makelink{param:simpleq-hardcore_xmin}{}
+    {\tt xmin} ({\tt Float64}, default: 0): minimum of the range of $|x|$ to be printed.
+
+    \item\makelink{param:simpleq-hardcore_xmax}{}
+    {\tt xmax} ({\tt Float64}, default: 100): maximum of the range of $|x|$ to be printed.
+
+    \item\makelink{param:simpleq-hardcore_nx}{}
+    {\tt nx} ({\tt Int64}, default: 100): number of points to print (linearly spaced).
+
+  \end{itemize}
+  \underline{Output} (one line for each value of $x$): [$|x|$] [$u(|x|)$]
+
+\end{itemize}
+
+\subsubsection{Description}
+\indent
+In order to carry out the computation of the solution of\-~(\ref{linearhc}) and compute the condensate fraction at the same time, we will consider the equation with an added parameter $\mu>0$:
+\begin{equation}
+  (-\Delta+4\epsilon)u=2e\rho u\ast u
+  ,\quad
+  \epsilon:=e+\frac\mu2
+  \label{hardcore_simple}
+\end{equation}
+for $|x|>1$.
+\bigskip
+
+\point{\bf Energy.}
+To compute the energy $e$ of this equation, with the extra parameter $\mu$, we consider the limit of the soft sphere potential $\lambda\mathds 1_{|x|<1}$ (see\-~(\ref{easyeq}) with $\beta_K=\beta_L=0$):
+\begin{equation}
+  (-\Delta+2\mu+4e)u(x)=s_\lambda(x)+2e\rho u\ast u
+  ,\quad
+  s_\lambda(x):=\lambda(1-u(x))\mathds 1_{|x|\leqslant 1}
+  ,\quad
+  \frac{2e}\rho=\int dx\ s_\lambda(x)
+  .
+\end{equation}
+Furthermore, since $\partial u$ need not be continuous at $|x|=1$, by integrating $-\Delta u$ over a thin spherical shell of radius 1, we find that, for $|x|\leqslant 1$,
+\begin{equation}
+  -\Delta u(x)=-\delta(|x|-1)\partial u|_{|x|\searrow1}
+\end{equation}
+so, formally,
+\begin{equation}
+  s_\infty(x)=\mathds 1_{|x|\leqslant 1}(2e(2-\rho u\ast u(x))+2\mu)-\delta(|x|-1)\partial u|_{|x|\searrow1}
+\end{equation}
+and
+\begin{equation}
+  \frac{2e}\rho=\int dx\ s_\infty(x)
+  =
+  2e\left(\frac{8\pi}3-\rho\int_{|x|<1}u\ast u(x)\right)+\frac{8\pi}3\mu-4\pi\partial u|_{|x|\searrow1}
+  .
+\end{equation}
+Therefore,
+\begin{equation}
+  e=\frac{2\pi\rho(\frac23\mu-\partial u|_{|x|\searrow1})}{1-\frac{8\pi}3\rho+\rho^2\int_{|x|<1}dx\ u\ast u(x)}
+  .
+  \label{emu}
+\end{equation}
+\bigskip
+
+\point{\bf Integral equation.}
+We turn the differential equation in\-~(\ref{hardcore_simple}) into an integral equation.
+Let
+\begin{equation}
+  w(|x|):=|x|u(x)
+\end{equation}
+we have, for $r>1$,
+\begin{equation}
+  (-\partial_r^2+4\epsilon)w(r)=2e\rho ru\ast u(r)
+  .
+\end{equation}
+Furthermore, the bounded solution of
+\begin{equation}
+  (-\partial_r^2+\alpha^2)w(r)=F(r)
+  ,\quad
+  w(1)=1
+\end{equation}
+is
+\begin{equation}
+  w(r)=e^{-\alpha(r-1)}+\int_0^\infty ds\ \frac{F(s+1)}{2\alpha}\left(e^{-\alpha|(r-1)-s|}-e^{-\alpha((r-1)+s)}\right)
+  \label{solw}
+\end{equation}
+so, for $r>1$,
+\begin{equation}
+  u(r)=
+  \frac 1re^{-2\sqrt\epsilon(r-1)}
+  +
+  \frac{\rho e}{2r\sqrt\epsilon}\int_0^\infty ds\ (s+1)(u\ast u(s+1))\left(e^{-2\sqrt\epsilon|(r-1)-s|}-e^{-2\sqrt\epsilon((r-1)+s)}\right)
+  .
+\end{equation}
+In order to compute the integral more easily, we split it:
+\begin{equation}
+  \begin{largearray}
+    u(r)=
+    \frac 1re^{-2\sqrt\epsilon(r-1)}
+    +
+    \frac{\rho e}{r\sqrt\epsilon}\int_0^{r-1} ds\ (s+1)(u\ast u(s+1))\sinh(2\sqrt\epsilon s)e^{-2\sqrt\epsilon(r-1)}
+    \\[0.5cm]\hfill
+    +
+    \frac{\rho e}{r\sqrt\epsilon}\sinh(2\sqrt\epsilon(r-1))\int_{r-1}^\infty ds\ (s+1)(u\ast u(s+1))e^{-2\sqrt\epsilon s}
+    .
+  \end{largearray}
+\end{equation}
+We change variables in the last integral:
+\begin{equation}
+  \begin{largearray}
+    u(r)=
+    \frac 1re^{-2\sqrt\epsilon(r-1)}
+    +
+    \frac{\rho e}{r\sqrt\epsilon}\int_0^{r-1} ds\ (s+1)(u\ast u(s+1))\sinh(2\sqrt\epsilon s)e^{-2\sqrt\epsilon(r-1)}
+    \\[0.5cm]\hfill
+    +
+    \frac{\rho e}{2r\sqrt\epsilon}\left(1-e^{-4\sqrt\epsilon(r-1)}\right)\int_0^\infty d\sigma\ (\sigma+r)(u\ast u(\sigma+r))e^{-2\sqrt\epsilon\sigma}
+    .
+  \end{largearray}
+  \label{uinteq}
+\end{equation}
+\bigskip
+
+\point{\bf The auto-convolution term.} We split
+\begin{equation}
+  u(r)=\mathds 1_{r>1}u_+(r-1)+\mathds 1_{r\leqslant 1}
+  .
+\end{equation}
+in terms of which
+\begin{equation}
+  u\ast u=
+  \mathds 1_{r\leqslant 1}\ast\mathds 1_{r\leqslant 1}
+  +2\mathds 1_{r\leqslant 1}\ast(u_+(r-1)\mathds 1_{r>1})
+  +(\mathds 1_{r>1}u_+(r-1))\ast(u_+(r-1)\mathds 1_{r>1})
+  \label{u*u_dcmp}
+\end{equation}
+In bipolar coordinates (see lemma\-~\ref{lemma:bipolar}),
+\begin{equation}
+  \mathds 1_{r\leqslant 1}\ast\mathds 1_{r\leqslant 1}(r)
+  =
+  \frac{2\pi}{r}\int_0^1 dt\ t\int_{|r-t|}^{r+t}ds\ s\mathds 1_{s\leqslant 1}
+  \label{1*1}
+\end{equation}
+and, if $r\geqslant 0$ and $t\geqslant 0$, then
+\begin{equation}
+  \int_{|r-t|}^{r+t}ds\ s\mathds 1_{s\leqslant 1}
+  =
+  \mathds 1_{r-1\leqslant t\leqslant r+1}\int_{|r-t|}^{\min(1,r+t)}ds\ s
+  =
+  \mathds 1_{r-1\leqslant t\leqslant r+1}
+  \left(
+    \mathds 1_{t\leqslant 1-r}
+    2rt
+    +
+    \mathds 1_{t>1-r}
+    \frac{1-(r-t)^2}2
+  \right)
+  .
+  \label{int1}
+\end{equation}
+Therefore, if $r\geqslant 1$
+\begin{equation}
+  \mathds 1_{r\leqslant 1}\ast\mathds 1_{r\leqslant 1}(r)
+  =
+  \mathds 1_{r\leqslant 2}
+  \frac{2\pi}{r}\int_{r-1}^1 dt\ t
+  \frac{1-(r-t)^2}2
+  =
+  \mathds 1_{r\leqslant 2}\frac\pi{12}(r-2)^2(r+4)
+  \label{u*u_1}
+\end{equation}
+and if $r<1$,
+\begin{equation}
+  \mathds 1_{r\leqslant 1}\ast\mathds 1_{r\leqslant 1}(r)
+  =
+  \frac{2\pi}r\int_0^1dt\ t
+  \left(
+    \mathds 1_{t\leqslant 1-r}
+    2rt
+    +
+    \mathds 1_{t>1-r}
+    \frac{1-(r-t)^2}2
+  \right)
+\end{equation}
+so
+\begin{equation}
+  \mathds 1_{r\leqslant 1}\ast\mathds 1_{r\leqslant 1}(r)
+  =\frac43\pi(1-r)^3
+  +\frac\pi{12}r(36-48r+17r^2)
+  =
+  \frac\pi{12}(r-2)^2(r+4)
+  .
+\end{equation}
+Thus, (\ref{u*u_1}) holds for all $r$.
+Furthermore,
+\begin{equation}
+  2\mathds 1_{r\leqslant 1}\ast(u_+(r-1)\mathds 1_{r>1})
+  =
+  \frac{4\pi}{r}\int_1^\infty dt\ tu_+(t-1)\int_{|r-t|}^{r+t}ds\ s\mathds 1_{s\leqslant 1}
+\end{equation}
+so, if $r\geqslant 0$ then, by\-~(\ref{int1}),
+\begin{equation}
+  2\mathds 1_{r\leqslant 1}\ast(u_+(r-1)\mathds 1_{r>1})(r)
+  =
+  \frac{2\pi}{r}\int_{\max(1,r-1)}^{r+1} dt\ tu_+(t-1)(1-(r-t)^2)
+  .
+  \label{u*u_cross}
+\end{equation}
+Finally, if $r\geqslant 0$ then
+\begin{equation}
+  (\mathds 1_{r>1}u_+(r-1))\ast(u_+(r-1)\mathds 1_{r>1})(r)
+  =
+  \frac{2\pi}{r}\int_1^\infty dt\ tu_+(t-1)\int_{\max{(1,|r-t|)}}^{r+t}ds\ su_+(s-1)
+  .
+  \label{u*u_u}
+\end{equation}
+Thus, by~\-(\ref{u*u_1}), (\ref{u*u_cross}) and\-~(\ref{u*u_u}), for $r\geqslant -1$,
+\begin{equation}
+  \begin{largearray}
+    u\ast u(1+r)=
+    \mathds 1_{r\leqslant 1}\frac\pi{12}(r-1)^2(r+5)
+    +
+    \frac{2\pi}{r+1}\int_{\max(0,r-1)}^{r+1} dt\ (t+1)u_+(t)(1-(r-t)^2)
+    \\[0.5cm]\hfill
+    +
+    \frac{2\pi}{r+1}\int_0^\infty dt\ (t+1)u_+(t)\int_{\max{(0,|r-t|-1)}}^{r+t+1}ds\ (s+1)u_+(s)
+    .
+  \end{largearray}
+\end{equation}
+We then compactify the integrals by changing variables to $\tau=\frac{1-t}{1+t}$ and $\sigma=\frac{1-s}{1+s}$:
+\begin{equation}
+  \begin{largearray}
+    u\ast u(1+r)=
+    \mathds 1_{r\leqslant 1}\frac\pi{12}(r-1)^2(r+5)
+    +
+    \frac{8\pi}{r+1}\int_{-\frac r{2+r}}^{\min(1,\frac{2-r}r)} d\tau\ \frac1{(1+\tau)^3}u_+({\textstyle\frac{1-\tau}{1+\tau}})(1-(r-{\textstyle\frac{1-\tau}{1+\tau}})^2)
+    \\[0.5cm]\hfill
+    +
+    \frac{32\pi}{r+1}\int_{-1}^1 d\tau\ \frac1{(1+\tau)^3}u_+({\textstyle\frac{1-\tau}{1+\tau}})\int_{\alpha_-(1+r,\tau)}^{\min(1,\beta_+(r,\tau))}d\sigma\ \frac1{(1+\sigma)^3}u_+({\textstyle\frac{1-\sigma}{1+\sigma}})
+    .
+  \end{largearray}
+  \label{u*u}
+\end{equation}
+with
+\begin{equation}
+  \alpha_-(1+r,\tau):=\frac{-r-\frac{1-\tau}{1+\tau}}{2+r+\frac{1-\tau}{1+\tau}}
+  ,\quad
+  \beta_+(r,\tau):=\frac{2-|r-\frac{1-\tau}{1+\tau}|}{|r-\frac{1-\tau}{1+\tau}|}
+\end{equation}
+(note that $\alpha_-$ is the same as defined for fulleq).
+Finally, note that, if $\beta_+<1$, that is, if $|r-\frac{1-\tau}{1+\tau}|>1$, then
+\begin{equation}
+  \beta_+(r,\tau)=\alpha_+({\textstyle |r-\frac{1-\tau}{1+\tau}|-1+\frac{1-\tau}{1+\tau}},\tau)
+\end{equation}
+where
+\begin{equation}
+  \alpha_+(r,\tau):=\frac{1-|r-\frac{1-\tau}{1+\tau}|}{1+|r-\frac{1-\tau}{1+\tau}|}
+\end{equation}
+is the same as is defined for \refname{sec:anyeq}{{\tt anyeq}} (\ref{anyeq_alpha}).
+\bigskip
+
+\point{\bf Chebyshev polynomial expansion.} We use the same interpolation as we used in \refname{sec:anyeq}{{\tt anyeq}}: (\ref{a_chebyshev_spline})
+\begin{equation}
+  \frac1{(1+\tau)^{\nu_u}}u_+({\textstyle\frac{1-\tau}{1+\tau}})=\sum_{l=0}^{J-1}\mathds 1_{\tau_l\leqslant\tau<\tau_{l+1}}\sum_{n=0}^\infty \mathbf F_{l,n}^{(\nu_u)}(u_+) T_n({\textstyle\frac{2\tau-(\tau_l+\tau_{l+1})}{\tau_{l+1}-\tau_l}})
+  \label{a_chebyshev_spline}
+\end{equation}
+($J$ is set by the parameter \refname{param:simpleq-hardcore_J}{{\tt J}})
+with
+\begin{equation}
+  \mathbf F_{l,n}^{(\nu_u)}(u_+):=
+  \frac{2-\delta_{n,0}}{\pi}
+  \int_0^{\pi} d\theta\ 
+  \frac{\cos(n\theta)}{(1+\frac{\tau_{l+1}-\tau_l}2\cos(\theta)+\frac{\tau_{l+1}+\tau_l}2)^{\nu_u}}
+  u_+({\textstyle\frac{2-(\tau_{l+1}-\tau_l)\cos(\theta)-(\tau_{l+1}+\tau_\alpha)}{2+(\tau_{l+1}-\tau_l)\cos(\theta)+(\tau_{l+1}+\tau_l)}})
+  \label{bfF}
+\end{equation}
+and we take $\nu_u$ to be the decay exponent of $u$, which we will assume is $\nu_u=4$.
+In particular, by\-~(\ref{u*u}),
+\begin{equation}
+  u\ast u(1+r)=\mathds 1_{r\leqslant 1}B^{(0)}(r)+\sum_l\sum_n\mathbf F_{l,n}^{(\nu_u)}(u_+)B^{(1)}_{l,n}(r)+\sum_{l,l'}\sum_{n,m}\mathbf F_{l,n}^{(\nu_u)}(u_+)\mathbf F_{l',m}^{(\nu_u)}(u_+)B^{(2)}_{l,n;l',m}(r)
+  \label{u*u}
+\end{equation}
+with
+\begin{equation}
+  B^{(0)}(r):=
+  \frac\pi{12}(r-1)^2(r+5)
+\end{equation}
+\begin{equation}
+  B^{(1)}_{l,n}(r):=
+  \mathds 1_{\tau_l<\frac{2-r}r}\mathds 1_{\tau_{l+1}>-\frac r{2+r}}
+  \frac{8\pi}{r+1}\int_{\max(\tau_l,-\frac r{2+r})}^{\min(\tau_{l+1},\frac{2-r}r)} d\tau\ \frac1{(1+\tau)^{3-\nu_u}}(1-(r-{\textstyle\frac{1-\tau}{1+\tau}})^2)T_n({\textstyle\frac{2\tau-(\tau_l+\tau_{l+1})}{\tau_{l+1}-\tau_l}})
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    B^{(2)}_{l,n;l',m}(r):=
+    \frac{32\pi}{r+1}\int_{\tau_l}^{\tau_{l+1}} d\tau\ \frac1{(1+\tau)^{3-\nu_u}}T_n({\textstyle\frac{2\tau-(\tau_l+\tau_{l+1})}{\tau_{l+1}-\tau_l}})
+    \mathds 1_{\tau_{l'}<\alpha_+(|r-\frac{1-\tau}{1+\tau}|-\frac{2\tau}{1+\tau},\tau)}\mathds 1_{\tau_{l'+1}>\alpha_-(1+r,\tau)}
+    \cdot\\[0.5cm]\hfill\cdot
+    \int_{\max(\tau_{l'},\alpha_-(1+r,\tau))}^{\min(\tau_{l'+1},\alpha_+(|r-\frac{1-\tau}{1+\tau}|-\frac{2\tau}{1+\tau},\tau))}d\sigma\ \frac1{(1+\sigma)^{3-\nu_u}}T_m({\textstyle\frac{2\sigma-(\tau_{l'}+\tau_{l'+1})}{\tau_{l'+1}-\tau_{l'}}})
+    .
+  \end{largearray}
+\end{equation}
+Thus, by\-~(\ref{uinteq}), for $r>0$,
+\begin{equation}
+  u_+(r,e,\epsilon)=
+  D^{(0)}(r,e,\epsilon)
+  +
+  \sum_l\sum_n\mathbf F_{l,n}^{(\nu_u)}(u_+)D^{(1)}_{l,n}(r,e,\epsilon)
+  +
+  \sum_{l,l'}\sum_{n,m}\mathbf F_{l,n}^{(\nu_u)}(u_+)\mathbf F_{l',m}^{(\nu_u)}(u_+)D^{(2)}_{l,n;l',m}(r,e,\epsilon)
+  \label{hardcore_u+}
+\end{equation}
+with
+\begin{equation}
+  \begin{largearray}
+    D^{(0)}(r,e,\epsilon):=
+    \frac 1{r+1}e^{-2\sqrt\epsilon r}
+    +
+    \frac{\rho e}{(r+1)\sqrt\epsilon}\int_0^{\min(1,r)} ds\ (s+1)B^{(0)}(s)\sinh(2\sqrt\epsilon s)e^{-2\sqrt\epsilon r}
+    \\[0.5cm]\hfill
+    +
+    \mathds 1_{r\leqslant 1}
+    \frac{\rho e}{2(r+1)\sqrt\epsilon}\left(1-e^{-4\sqrt\epsilon r}\right)\int_0^{1-r} d\sigma\ (\sigma+r+1)B^{(0)}(\sigma+r)e^{-2\sqrt\epsilon\sigma}
+  \end{largearray}
+  \label{D0}
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    D^{(1)}_{l,n}(r,e,\epsilon):=
+    \frac{\rho e}{(r+1)\sqrt\epsilon}\int_0^r ds\ (s+1)B^{(1)}_{l,n}(s)\sinh(2\sqrt\epsilon s)e^{-2\sqrt\epsilon r}
+    \\[0.5cm]\hfill
+    +
+    \frac{\rho e}{2(r+1)\sqrt\epsilon}\left(1-e^{-4\sqrt\epsilon r}\right)\int_0^\infty d\sigma\ (\sigma+r+1)B_{l,n}^{(1)}(\sigma+r)e^{-2\sqrt\epsilon\sigma}
+    .
+  \end{largearray}
+  \label{D1}
+\end{equation}
+and
+\begin{equation}
+  \begin{largearray}
+    D^{(2)}_{l,n;l',m}(r,e,\epsilon):=
+    \frac{\rho e}{(r+1)\sqrt\epsilon}\int_0^r ds\ (s+1)B^{(2)}_{l,n;l',m}(s)\sinh(2\sqrt\epsilon s)e^{-2\sqrt\epsilon r}
+    \\[0.5cm]\hfill
+    +
+    \frac{\rho e}{2(r+1)\sqrt\epsilon}\left(1-e^{-4\sqrt\epsilon r}\right)\int_0^\infty d\sigma\ (\sigma+r+1)B_{l,n;l',m}^{(2)}(\sigma+r)e^{-2\sqrt\epsilon\sigma}
+    .
+  \end{largearray}
+  \label{D2}
+\end{equation}
+\bigskip
+
+\point{\bf Energy.} We now compute the approximation for the energy, using\-~(\ref{hardcore_energy}).
+\bigskip
+
+\subpoint We start with $\partial u|_{|x|\searrow 1}$.
+By\-~(\ref{uinteq}),
+\begin{equation}
+  \partial u|_{|x|\searrow 1}
+  =
+  -(1+2\sqrt\epsilon)
+  +
+  2\rho e\int_0^\infty d\sigma\ (\sigma+1)(u\ast u(\sigma+1))e^{-2\sqrt\epsilon\sigma}
+\end{equation}
+\bigskip
+which, by\-~(\ref{u*u}), becomes
+\begin{equation}
+  \begin{largearray}
+    \partial u|_{|x|\searrow 1}
+    =
+    -(1+2\sqrt\epsilon)
+    +
+    \gamma^{(0)}(e,\epsilon)
+    +
+    \sum_l\sum_n\mathbf F_{l,n}^{(\nu_u)}(u_+)\gamma^{(1)}_{l,n}(e,\epsilon)
+    +\\\hfill+
+    \sum_{l,l'}\sum_{n,m}\mathbf F_{l,n}^{(\nu_u)}(u_+)\mathbf F_{l',m}^{(\nu_u)}(u_+)\gamma^{(2)}_{l,n;l',m}(e,\epsilon)
+  \end{largearray}
+\end{equation}
+with
+\begin{equation}
+  \gamma^{(0)}(e,\epsilon):=
+  2\rho e\int_0^1 d\sigma\ (\sigma+1)B^{(0)}(\sigma)e^{-2\sqrt\epsilon\sigma}
+\end{equation}
+\begin{equation}
+  \gamma^{(1)}_{l,n}(e,\epsilon):=
+  2\rho e\int_0^\infty d\sigma\ (\sigma+1)B_{l,n}^{(1)}(\sigma)e^{-2\sqrt\epsilon\sigma}
+\end{equation}
+and
+\begin{equation}
+  \gamma^{(2)}_{l,n;l',m}(e,\epsilon):=
+  2\rho e\int_0^\infty d\sigma\ (\sigma+1)B_{l,n;l',m}^{(2)}(\sigma)e^{-2\sqrt\epsilon\sigma}
+  .
+\end{equation}
+
+\subpoint Let us now turn to $\int_{|x|<1}dx\ u\ast u(x)$.
+We have
+\begin{equation}
+  \int_{|x|<1}dx\ u\ast u(x)
+  =
+  4\pi\int_0^1 dr\ r^2 u\ast u(r)
+\end{equation}
+so, by\-~(\ref{u*u}),
+\begin{equation}
+  \int_{|x|<1}dx\ u\ast u(x)
+  =
+  \bar\gamma^{(0)}(r)
+  +
+  \sum_l\sum_n\mathbf F_{l,n}^{(\nu_u)}(u_+)\bar\gamma^{(1)}_{l,n}(r)
+  +
+  \sum_{l,l'}\sum_{n,m}\mathbf F_{l,n}^{(\nu_u)}(u_+)\mathbf F_{l',m}^{(\nu_u)}(u_+)\bar\gamma^{(2)}_{l,n;l',m}(r)
+  \label{intusu}
+\end{equation}
+with
+\begin{equation}
+  \bar\gamma^{(0)}:=
+  4\pi\int_0^1 d\sigma\ \sigma^2B^{(0)}(\sigma-1)
+\end{equation}
+\begin{equation}
+  \bar\gamma^{(1)}_{l,n}:=
+  4\pi\int_0^1 d\sigma\ \sigma^2B_{l,n}^{(1)}(\sigma-1)
+\end{equation}
+and
+\begin{equation}
+  \bar\gamma^{(2)}_{l,n;l',m}:=
+  4\pi\int_0^1 d\sigma\ \sigma^2B_{l,n;l',m}^{(2)}(\sigma-1)
+  .
+\end{equation}
+\bigskip
+
+\subpoint Thus,
+\begin{equation}
+  \begin{largearray}
+    e=2\pi\rho
+    \cdot\\\hfill\cdot
+    \frac
+    {
+      C(\epsilon)
+      -
+      \gamma^{(0)}(e,\epsilon)
+      -
+      \sum_l\sum_n\mathbf F_{l,n}^{(\nu_u)}(u_+)\gamma^{(1)}_{l,n}(e,\epsilon)
+      -
+      \sum_{l,l'}\sum_{n,m}\mathbf F_{l,n}^{(\nu_u)}(u_+)\mathbf F_{l',m}^{(\nu_u)}(u_+)\gamma^{(2)}_{l,n;l',m}(e,\epsilon)
+    }
+    {1-\frac83\pi\rho+\rho^2\left(
+      \bar\gamma^{(0)}
+      +
+      \sum_l\sum_n\mathbf F_{l,n}^{(\nu_u)}(u_+)\bar\gamma^{(1)}_{l,n}
+      +
+      \sum_{l,l'}\sum_{n,m}\mathbf F_{l,n}^{(\nu_u)}(u_+)\mathbf F_{l',m}^{(\nu_u)}(u_+)\bar\gamma^{(2)}_{l,n;l',m}
+    \right)}
+    .
+  \end{largearray}
+\end{equation}
+\begin{equation}
+  C(\epsilon):=\frac23\mu+(1+2\sqrt\epsilon)
+  =\frac43(\epsilon-e)+1+2\sqrt\epsilon
+\end{equation}
+\bigskip
+
+\point{\bf Newton algorithm.}
+In this paragraph, we set $\epsilon=e$, that is, $\mu=0$.
+\bigskip
+
+\subpoint
+As we did for \refname{sec:anyeq}{{\tt anyeq}}, we discretize the integral in\-~(\ref{bfF}) by using a Gauss-Legendre quadrature, and truncate the sum over Chebyshev polynomials to order \refname{param:simpleq-hardcore_P}{{\tt P}}.
+We then reduce the computation to a finite system of equations, whose variables are
+\begin{equation}
+  e
+  ,\quad
+  \mathfrak u_{l,j}:=u_+(r_{l,j})
+  ,\quad
+  r_{l,j}:=
+  \frac{2+(\tau_{l+1}-\tau_l)\sin(\frac{\pi x_j}2)-(\tau_{l+1}+\tau_\alpha)}{2-(\tau_{l+1}-\tau_l)\sin(\frac{\pi x_j}2)+(\tau_{l+1}+\tau_l)}
+\end{equation}
+where $x_j$ are the abcissa for Gauss-Legendre quadratures (see\-~(\ref{klj})).
+In other words, we define a vector $\mathbb U$ in dimension \refname{param:simpleq-hardcore_N}{{\tt N}}$\times$\refname{param:simpleq-hardcore_J}{{\tt J}}$+1$, where the first \refname{param:simpleq-hardcore_N}{{\tt N}}$\times$\refname{param:simpleq-hardcore_J}{{\tt J}} terms are $\mathfrak u_{l,j}$ and the last component is $e$.
+We then write\-~(\ref{hardcore_u+}) as, for $l\in\{0,\cdots,J-1\}$, $j\in\{1,\cdots,N\}$,
+\begin{equation}
+  \Xi_{lN+j}(\mathbb U)=0
+  ,\quad
+  \Xi_{NJ+1}(\mathbb U)=0
+\end{equation}
+(note that $\Xi_{lN+j}$ corresponds to the pair $(l,j)$)
+with
+\begin{equation}
+  \begin{largearray}
+    \Xi_{lN+j}(\mathbb U):=
+    -\mathfrak u_{l,j}
+    +
+    \mathbb D^{(0)}(r_{l,j},e)
+    +
+    \sum_{l'}\sum_n\mathfrak F_{l',n}^{(\nu_u)}(\mathfrak u)\mathbb D^{(1)}_{l',n}(r_{l,j},e)
+    +\\\hfill+
+    \sum_{l',l''}\sum_{n,m}\mathfrak F_{l',n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathfrak u)\mathbb D^{(2)}_{l',n;l'',m}(r_{l,j},e)
+  \end{largearray}
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \Xi_{NJ+1}(\mathbb U):=-e
+    +\\\hfill+
+    2\pi\rho
+    \frac
+    {
+      C(e)
+      -
+      \mathfrak g^{(0)}(e)
+      -
+      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak g^{(1)}_{l,n}(e)
+      -
+      \sum_{l,l'}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l',m}^{(\nu_u)}(\mathfrak u)\mathfrak g^{(2)}_{l,n;l',m}(e)
+    }
+    {1-\frac83\pi\rho+\rho^2\left(
+      \bar{\mathfrak g}^{(0)}
+      +
+      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(1)}_{l,n}
+      +
+      \sum_{l,l'}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l',m}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(2)}_{l,n;l',m}
+    \right)}
+    .
+  \end{largearray}
+\end{equation}
+where $\mathfrak F$ is defined in\-~(\ref{frakF}), and $\mathbb D^{(i)}$, $\mathfrak g^{(i)}$, $\bar{\mathfrak g}^{(i)}$ are defined like $D^{(i)}$, $\gamma^{(i)}$ and $\bar\gamma^{(i)}$ except that the integrals over bounded intervals are approximated using Gauss-Legendre quadratures, and the integrals from $0$ to $\infty$ are approximated using Gauss-Laguerre quadratures.
+Gauss-Legendre and Gauss-Laguerre quadratures and their errors are discussed in appendix\-~\ref{appGL}.
+The orders of the quadratures are given by the variable \refname{param:simpleq-hardcore_N}{{\tt N}}.
+\bigskip
+
+\subpoint
+We then solve $\Xi=0$ using the Newton algorithm, that is, we define a sequence of $\mathbb U$'s:
+\begin{equation}
+  \mathbb U^{(n+1)}=\mathbb U^{(n)}-(D\Xi(\mathbb U^{(n)}))^{-1}\Xi(\mathbb U^{(n)})
+\end{equation}
+where $D\Xi$ is the Jacobian of $\Xi$:
+\begin{equation}
+  (D\Xi)_{\alpha,\beta}:=\frac{\partial\Xi_\alpha}{\partial\mathbb U_\beta}
+  .
+\end{equation}
+We initialize the algorithm with
+\begin{equation}
+  \mathbb U^{(0)}_{lN+j}=\frac1{(1+r_{l,j}^2)^2}
+  ,\quad
+  \mathbb U^{(0)}_{JN+1}=\pi\rho
+  .
+\end{equation}
+\bigskip
+
+\subpoint
+We are left with computing the Jacobian of $\Xi$:
+\begin{equation}
+  \begin{largearray}
+    \frac{\partial\Xi_{lN+j}(\mathbb U)}{\partial\mathfrak u_{l',i}}=
+    -\delta_{l,l'}\delta_{j,i}
+    +
+    \sum_{l''}\sum_n\mathfrak F_{l'',n}^{(\nu_u)}(\mathds 1_{l',i})\mathbb D^{(1)}_{l'',n}(r_{l,j},e)
+    +\\\hfill+
+    2\sum_{l'',l'''}\sum_{n,m}\mathfrak F_{l'',n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l''',m}^{(\nu_u)}(\mathds 1_{l',j})\mathbb D^{(2)}_{l'',n;l''',m}(r_{l,j},e)
+  \end{largearray}
+\end{equation}
+where $\mathds 1_{l',i}$ is a vector whose only non-vanishing entry is the $(l',i)$-th, which is 1,
+\begin{equation}
+  \begin{largearray}
+    \frac{\partial\Xi_{JN+1}(\mathbb U)}{\partial\mathfrak u_{l',i}}
+    =\\\hfill=
+    2\pi\rho
+    \frac
+    {
+      -
+      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathds 1_{l',i})\mathfrak g^{(1)}_{l,n}(e)
+      -
+      2\sum_{l,l''}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathds 1_{l',i})\mathfrak g^{(2)}_{l,n;l'',m}(e)
+    }
+    {1-\frac83\pi\rho+\rho^2\left(
+      \bar{\mathfrak g}^{(0)}
+      +
+      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(1)}_{l,n}
+      +
+      \sum_{l,l''}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(2)}_{l,n;l'',m}
+    \right)}
+    \\[1cm]\hfill
+    -(\Xi_{NJ+1}(\mathbb U)+e)
+    \frac
+    {\rho^2\left(
+      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathds 1_{l'.i})\bar{\mathfrak g}^{(1)}_{l,n}
+      +
+      2\sum_{l,l''}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathds 1_{l',i})\bar{\mathfrak g}^{(2)}_{l,n;l'',m}
+    \right)}{1-\frac83\pi\rho+\rho^2\left(
+      \bar{\mathfrak g}^{(0)}
+      +
+      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(1)}_{l,n}
+      +
+      \sum_{l,l''}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(2)}_{l,n;l'',m}
+    \right)}
+    .
+  \end{largearray}
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \frac{\partial\Xi_{lN+j}(\mathbb U)}{\partial e}=
+    \partial_e\mathbb D^{(0)}(r_{l,j},e)
+    +
+    \sum_{l'}\sum_n\mathfrak F_{l',n}^{(\nu_u)}(\mathfrak u)\partial_e\mathbb D^{(1)}_{l',n}(r_{l,j},e)
+    +\\\hfill+
+    \sum_{l',l''}\sum_{n,m}\mathfrak F_{l',n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathfrak u)\partial_e\mathbb D^{(2)}_{l',n;l'',m}(r_{l,j},e)
+  \end{largearray}
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \frac{\partial\Xi_{NJ+1}(\mathbb U)}{\partial e}
+    =
+    -1
+    +\\\hfill+
+    2\pi\rho
+    \frac
+    {
+      \partial_e C(e)
+      -
+      \partial_e\mathfrak g^{(0)}(e)
+      -
+      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\partial_e\mathfrak g^{(1)}_{l,n}(e)
+      -
+      \sum_{l,l'}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l',m}^{(\nu_u)}(\mathfrak u)\partial_e\mathfrak g^{(2)}_{l,n;l',m}(e)
+    }
+    {1-\frac83\pi\rho+\rho^2\left(
+      \bar{\mathfrak g}^{(0)}
+      +
+      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(1)}_{l,n}
+      +
+      \sum_{l,l'}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l',m}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(2)}_{l,n;l',m}
+    \right)}
+    .
+  \end{largearray}
+\end{equation}
+To compute $\partial_e\mathbb D^{(i)}$ and $\partial_e\mathfrak g^{(i)}$, we use $\partial_e=\frac1{2\sqrt e}\frac\partial{\partial\sqrt e}$ and
+\begin{equation}
+  \frac{\partial C(e)}{\partial\sqrt e}=2
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \frac{\partial D^{(0)}(r,e)}{\partial\sqrt e}
+    =
+    -\frac{2r}{r+1}e^{-2\sqrt er}
+    \\[0.5cm]\hfill
+    +
+    \frac{\rho}{r+1}\int_0^{\min(1,r)} ds\ (s+1)B^{(0)}(s)e^{-2\sqrt er}\left(
+      (1-2\sqrt er)\sinh(2\sqrt e s)
+      +2\sqrt es\cosh(2\sqrt e s)
+    \right)
+    \\[0.5cm]
+    +
+    \mathds 1_{r\leqslant 1}
+    \frac{\rho}{2(r+1)}\int_0^{1-r} d\sigma\ (\sigma+r+1)B^{(0)}(\sigma+r)
+    \cdot\\[0.5cm]\hfill\cdot
+    \left(
+      (1-2\sqrt e\sigma)\left(1-e^{-4\sqrt er}\right)e^{-2\sqrt e\sigma}
+      +
+      4\sqrt e re^{-2\sqrt e(2r+\sigma)}
+    \right)
+  \end{largearray}
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \frac{\partial D_{l,n}^{(1)}(r,e)}{\partial\sqrt e}
+    =
+    \frac{\rho}{r+1}\int_0^{r} ds\ (s+1)B_{l,n}^{(1)}(s)e^{-2\sqrt er}\left(
+      (1-2\sqrt er)\sinh(2\sqrt e s)
+      +2\sqrt es\cosh(2\sqrt e s)
+    \right)
+    \\[0.5cm]
+    +
+    \frac{\rho}{2(r+1)}\int_0^\infty d\sigma\ (\sigma+r+1)B_{l,n}^{(1)}(\sigma+r)
+    \cdot\\[0.5cm]\hfill\cdot
+    \left(
+      (1-2\sqrt e\sigma)\left(1-e^{-4\sqrt er}\right)e^{-2\sqrt e\sigma}
+      +
+      4\sqrt e re^{-2\sqrt e(2r+\sigma)}
+    \right)
+  \end{largearray}
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \frac{\partial D_{l,n;l',m}^{(2)}(r,e)}{\partial\sqrt e}
+    \\[0.5cm]
+    =
+    \frac{\rho}{r+1}\int_0^{r} ds\ (s+1)B_{l,n;l',m}^{(2)}(s)e^{-2\sqrt er}\left(
+      (1-2\sqrt er)\sinh(2\sqrt e s)
+      +2\sqrt es\cosh(2\sqrt e s)
+    \right)
+    \\[0.5cm]
+    +
+    \frac{\rho}{2(r+1)}\int_0^\infty d\sigma\ (\sigma+r+1)B_{l,n;l',m}^{(2)}(\sigma+r)
+    \cdot\\[0.5cm]\hfill\cdot
+    \left(
+      (1-2\sqrt e\sigma)\left(1-e^{-4\sqrt er}\right)e^{-2\sqrt e\sigma}
+      +
+      4\sqrt e re^{-2\sqrt e(2r+\sigma)}
+    \right)
+    .
+  \end{largearray}
+\end{equation}
+Furthermore,
+\begin{equation}
+  \frac{\partial\gamma^{(0)}(e)}{\partial\sqrt e}=
+  4\rho\sqrt e\int_0^1 d\sigma\ (\sigma+1)B^{(0)}(\sigma)(1-\sqrt e\sigma)e^{-2\sqrt e\sigma}
+\end{equation}
+\begin{equation}
+  \frac{\partial\gamma^{(1)}_{l,n}(e)}{\partial\sqrt e}=
+  4\rho e\int_0^\infty d\sigma\ (\sigma+1)B_{l,n}^{(1)}(\sigma)(1-\sqrt e\sigma)e^{-2\sqrt e\sigma}
+\end{equation}
+and
+\begin{equation}
+  \frac{\partial\gamma^{(2)}_{l,n;l',m}(e)}{\partial\sqrt e}=
+  4\rho e\int_0^\infty d\sigma\ (\sigma+1)B_{l,n;l',m}^{(2)}(\sigma)(1-\sqrt e\sigma)e^{-2\sqrt e\sigma}
+  .
+\end{equation}
+Finally, to get from $D$ to $\mathbb D$ and $\gamma$ to $\mathfrak g$, we approximate the integrate using Gauss-Legendre and Gauss-Laguerre quadratures (see appendix\-~\ref{appGL}), as described above.
+\bigskip
+
+\subpoint
+We iterate the Newton algorithm until the Newton relative error $\epsilon$ becomes smaller than the \refname{param:easyeq_tolerance}{{\tt tolerance}} parameter.
+The Newton error is defined as
+\begin{equation}
+  \epsilon=\frac{\|\mathbb U^{(n+1)}-\mathbb U^{(n)}\|_2}{\|\mathbb U^{(n)}\|_2}
+\end{equation}
+where $\|\cdot\|_2$ is the $l_2$ norm.
+The energy thus obtained is
+\begin{equation}
+  e=\mathbb U_{JN+1}
+  .
+\end{equation}
+\bigskip
+
+\point{\bf Condensate fraction.}
+To compute the condensate fraction, we use the parameter $\mu$ in\-~(\ref{hardcore_simple}).
+The uncondensed fraction is
+\begin{equation}
+  \eta=\partial_\mu e|_{\mu=0}
+  .
+\end{equation}
+To compute $\partial_\mu e$, we use
+\begin{equation}
+  \Xi(\mathbb U)=0
+\end{equation}
+which we differentiate with respect to $\mu$:
+\begin{equation}
+  \partial_\mu\mathbb U=-(D\Xi)^{-1}\partial_\mu\Xi
+  .
+\end{equation}
+We are left with computing $\partial_\mu\Xi$:
+\begin{equation}
+  \begin{largearray}
+    \frac{\partial\Xi_{lN+j}(\mathbb U)}{\partial\mu}=
+    \partial_\mu\mathbb D^{(0)}(r_{l,j},e,\epsilon)
+    +
+    \sum_{l'}\sum_n\mathfrak F_{l',n}^{(\nu_u)}(\mathfrak u)\partial_\mu\mathbb D^{(1)}_{l',n}(r_{l,j},e,\epsilon)
+    +\\\hfill+
+    \sum_{l',l''}\sum_{n,m}\mathfrak F_{l',n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathfrak u)\partial_\mu\mathbb D^{(2)}_{l',n;l'',m}(r_{l,j},e,\epsilon)
+  \end{largearray}
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \frac{\partial\Xi_{NJ+1}(\mathbb U)}{\partial \mu}
+    =\\\hfill=
+    2\pi\rho
+    \frac
+    {
+      \partial_\mu C(\epsilon)
+      -
+      \partial_\mu\mathfrak g^{(0)}(e,\epsilon)
+      -
+      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\partial_\mu\mathfrak g^{(1)}_{l,n}(e,\epsilon)
+      -
+      \sum_{l,l'}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l',m}^{(\nu_u)}(\mathfrak u)\partial_\mu\mathfrak g^{(2)}_{l,n;l',m}(e,\epsilon)
+    }
+    {1-\frac83\pi\rho+\rho^2\left(
+      \bar{\mathfrak g}^{(0)}
+      +
+      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(1)}_{l,n}
+      +
+      \sum_{l,l'}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l',m}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(2)}_{l,n;l',m}
+    \right)}
+    .
+  \end{largearray}
+\end{equation}
+We then use $\partial_\mu=\frac1{4\sqrt\epsilon}\partial_{\sqrt\epsilon}$ and
+\begin{equation}
+  \left.\frac{\partial C(\epsilon)}{\partial\sqrt\epsilon}\right|_{\mu=0}
+  =
+  \frac83\sqrt\epsilon+2
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \left.\frac{\partial D^{(0)}(r)}{\partial\sqrt\epsilon}\right|_{\mu=0}=
+    -\frac{2r}{r+1}e^{-2\sqrt{e}r}
+    \\[0.5cm]\hfill
+    +
+    \frac{\rho}{(r+1)}\int_0^{\min(1,r)} ds\ (s+1)B^{(0)}(s)e^{-2\sqrt{e}r}((-1-2\sqrt er)\sinh(2\sqrt{e} s)+2\sqrt es\cosh(2\sqrt{e} s))
+    \\[0.5cm]
+    +
+    \mathds 1_{r\leqslant 1}
+    \frac{\rho}{2(r+1)}\int_0^{1-r} d\sigma\ (\sigma+r+1)B^{(0)}(\sigma+r)
+    \cdot\\[0.5cm]\hfill\cdot
+    \left((-1-2\sqrt e\sigma)\left(1-e^{-4\sqrt er}\right)e^{-2\sqrt{e}\sigma}+4\sqrt ere^{-2\sqrt{e}(2r+\sigma)}\right)
+  \end{largearray}
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \left.\frac{\partial D_{l,n}^{(1)}(r)}{\partial\sqrt\epsilon}\right|_{\mu=0}
+    =
+    \frac{\rho}{r+1}\int_0^{r} ds\ (s+1)B_{l,n}^{(1)}(s)e^{-2\sqrt er}\left(
+      (-1-2\sqrt er)\sinh(2\sqrt e s)
+      +2\sqrt es\cosh(2\sqrt e s)
+    \right)
+    \\[0.5cm]
+    +
+    \frac{\rho}{2(r+1)}\int_0^\infty d\sigma\ (\sigma+r+1)B_{l,n}^{(1)}(\sigma+r)
+    \cdot\\[0.5cm]\hfill\cdot
+    \left(
+      (-1-2\sqrt e\sigma)\left(1-e^{-4\sqrt er}\right)e^{-2\sqrt e\sigma}
+      +
+      4\sqrt e re^{-2\sqrt e(2r+\sigma)}
+    \right)
+  \end{largearray}
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    \left.\frac{\partial D_{l,n;l',m}^{(2)}(r)}{\partial\sqrt\epsilon}\right|_{\mu=0}
+    \\[0.5cm]
+    =
+    \frac{\rho}{r+1}\int_0^{r} ds\ (s+1)B_{l,n;l',m}^{(2)}(s)e^{-2\sqrt er}\left(
+      (-1-2\sqrt er)\sinh(2\sqrt e s)
+      +2\sqrt es\cosh(2\sqrt e s)
+    \right)
+    \\[0.5cm]
+    +
+    \frac{\rho}{2(r+1)}\int_0^\infty d\sigma\ (\sigma+r+1)B_{l,n;l',m}^{(2)}(\sigma+r)
+    \cdot\\[0.5cm]\hfill\cdot
+    \left(
+      (-1-2\sqrt e\sigma)\left(1-e^{-4\sqrt er}\right)e^{-2\sqrt e\sigma}
+      +
+      4\sqrt e re^{-2\sqrt e(2r+\sigma)}
+    \right)
+    .
+  \end{largearray}
+\end{equation}
+\begin{equation}
+  \left.\frac{\partial \gamma^{(0)}}{\partial\sqrt\epsilon}\right|_{\mu=0}
+  =
+  -4\rho e\int_0^1 d\sigma\ \sigma(\sigma+1)B^{(0)}(\sigma)e^{-2\sqrt e\sigma}
+\end{equation}
+\begin{equation}
+  \left.\frac{\partial\gamma^{(1)}_{l,n}}{\partial\sqrt\epsilon}\right|_{\mu=0}
+  =
+  -4\rho e\int_0^\infty d\sigma\ (\sigma+1)\sigma B_{l,n}^{(1)}(\sigma)e^{-2\sqrt e\sigma}
+\end{equation}
+and
+\begin{equation}
+  \left.\frac{\partial\gamma^{(2)}_{l,n;l',m}}{\partial\sqrt\epsilon}\right|_{\mu=0}
+  =
+  -4\rho e\int_0^\infty d\sigma\ (\sigma+1)\sigma B_{l,n;l',m}^{(2)}(\sigma)e^{-2\sqrt e\sigma}
+  .
+\end{equation}
+
+\subsection{\tt simpleq-iteration}\label{sec:simpleq-iteration}
+\indent
+This method is used to solve the Simple equation using the iteration described in\-~\cite{CJL20}.
+The Simple equation is
+\begin{equation}
+  -\Delta u
+  =
+  S-4e u+2e\rho u\ast u
+  \label{simpleq_iteration}
+\end{equation}
+\begin{equation}
+  S:=(1-u)v
+  ,\quad
+  \rho:=\frac{2e}{\int dx\ (1-u(|x|))v(|x|)}
+  .
+  \label{simpleq_Se}
+\end{equation}
+for a soft potential $v$ at fixed energy $e>0$.
+The iteration is defined as
+\begin{equation}
+  -\Delta u_n=S_n-4eu_n+2e\rho_{n-1}u_{n-1}\ast u_{n-1}
+  ,\quad
+  u_0=0
+  ,\quad
+  \rho_{n-1}=
+  \frac{2e}{\int dx\ (1-u_{n-1}(|x|))v(|x|)}
+  .
+  \label{iteration}
+\end{equation}
+\bigskip
+
+\subsubsection{Usage}
+\indent
+Unless otherwise noted, this method takes the following parameters (specified via the {\tt [-p params]} flag, as a `{\tt;}' separated list of entries).
+\begin{itemize}
+  \item\makelink{param:simpleq-iteration_e}{}
+  {\tt e} ({\tt Float64}, default: $10^{-4}$): energy $e$.
+
+  \item\makelink{param:simpleq-iteration_maxiter}{}
+  {\tt maxiter} ({\tt Int64}, default: 21): maximal number of iterations.
+
+  \item\makelink{param:simpleq-iteration_order}{}
+  {\tt order} ({\tt Int64}, default: 100): order used for all Gauss quadratures (denoted by $N$ below).
+
+\end{itemize}
+\bigskip
+
+The available {\tt commands} are the following.
+
+\begin{itemize}
+  \item {\tt rho\_e}\makelink{command:simpleq-iteration_rho_e}{}:
+  compute the density $\rho$ as a function of $e$.\par
+  \underline{Disabled parameters}: \refname{param:simpleq-iteration_e}{{\tt e}}.\par
+  \underline{Extra parameters}:
+  \begin{itemize}
+    \item\makelink{param:simpleq-iteration_minle}{}
+    {\tt minle} ({\tt Float64}, default: $10^{-6}$): minimal value for $\log_{10}e$.
+
+    \item\makelink{param:simpleq-iteration_maxle}{}
+    {\tt maxle} ({\tt Float64}, default: $10^{2}$): maximal value for $\log_{10}e$.
+
+    \item\makelink{param:simpleq-iteration_nle}{}
+    {\tt nle} ({\tt Int64}, default: $100$): number of values for $e$ (spaced logarithmically).
+
+    \item\makelink{param:simpleq-iteration_es}{}
+    {\tt es} ({\tt Array\{Float64\}}, default: $(10^{{\tt minle}+\frac{{\tt maxle}-{\tt minle}}{{\tt nle}}n})_n$: list of values for $e$, specified as a `{\tt,}' separated list.
+  \end{itemize}
+  \underline{Output} (one line for each value of $e$): [$e$] [$\rho$].
+
+  \item\makelink{command:simpleq-iteration_ux}{}
+  {\tt ux}: compute $u$ as a function of $|x|$.\par
+  \underline{Extra parameters}:
+  \begin{itemize}
+    \item\makelink{param:simpleq-iteration_xmin}{}
+    {\tt xmin} ({\tt Float64}, default: 0): minimum of the range of $|x|$ to be printed.
+
+    \item\makelink{param:simpleq-iteration_xmax}{}
+    {\tt xmax} ({\tt Float64}, default: 100): maximum of the range of $|x|$ to be printed.
+
+    \item\makelink{param:simpleq-iteration_nx}{}
+    {\tt nx} ({\tt Int64}, default: 100): number of points to print (linearly spaced).
+  \end{itemize}
+  \underline{Output} (one line for each value of $x$): [$|x|$] [$u_1(|x|)$] [$u_2(|x|)$] [$u_3(|x|)$] $\cdots$
+
+\end{itemize}
+
+\subsubsection{Description}
+\indent
+In Fourier space
+\begin{equation}
+  \hat u_n(|k|):=\int dx\ e^{ikx}u_n(|x|)
+\end{equation}
+(\ref{iteration}) becomes
+\begin{equation}
+  (k^2+4e)\hat u_n(|k|)=\hat S_n(|k|)+2e\rho_{n-1}\hat u_{n-1}(k)^2
+  ,\quad
+  \rho_n:=\frac{2e}{\hat S_n(0)}
+  \label{iteration_uk}
+\end{equation}
+with
+\begin{equation}
+  \hat S_n(|k|)
+  =\hat v(|k|)-\frac1{(2\pi)^3}\int dp\ \hat u_n(|p|)\hat v(|k-p|)
+  .
+\end{equation}
+We write $\hat S_n$ in bipolar coordinates (see lemma\-~\ref{lemma:bipolar}):
+\begin{equation}
+  \hat S_n(|k|)=\hat v(|k|)-\frac1{(2\pi)^3}\int_0^\infty dt\ t\hat u_n(t)\Eta(|k|,t)
+\end{equation}
+with
+\begin{equation}
+  \Eta(y,t):=\frac{2\pi}y\int_{|y-t|}^{y+t}ds\ s\hat v(s)
+\end{equation}
+(note that this agrees with the function\-~(\ref{easyeq_Eta}) defined for \refname{sec:easyeq}{{\tt easyeq}}).
+We also change variables to
+\begin{equation}
+  y=\frac1{t+1}
+  ,\quad
+  t=\frac{1-y}y
+\end{equation}
+\begin{equation}
+  \hat S_n(|k|)=\hat v(|k|)-\frac1{(2\pi)^3}\int_0^1 dy\ \frac{(1-y)\hat u_n(\frac{1-y}y)\Eta(|k|,\frac{1-y}y)}{y^3}
+  .
+\end{equation}
+We approximate this integral using a Gauss-Legendre quadrature (see appendix\-~\ref{appGL}) and discretize Fourier space:
+\begin{equation}
+  k_i:=\frac{1-x_i}{1+x_i}
+  ,\quad
+  y_i:=\frac{x_i+1}2
+\end{equation}
+where $x_i$ are the Gauss-Legendre abscissa, and
+\begin{equation}
+  \hat S_n(k_i)\approx\hat v(k_i)-\frac1{2(2\pi)^3}\sum_{j=1}^Nw_j\frac{(1-y_j)\hat u_n(\frac{1-y_j}{y_j})\Eta(k_i,\frac{1-y_j}{y_j})}{y_j^3}
+\end{equation}
+so if
+\begin{equation}
+  \mathbb U_i(n):=\hat u_n(k_i)
+  ,\quad
+  \mathbf v_i:=\hat v(k_i)
+\end{equation}
+we have
+\begin{equation}
+  \sum_{j=1}^NA_{i,j}\mathbb U_j(n)=b_i^{(n)}
+\end{equation}
+with
+\begin{equation}
+  A_{i,j}:=(k_i^2+4e)\delta_{i,j}+\frac{w_j(1-y_j)\Eta(k_i,\frac{1-y_j}{y_j})}{2(2\pi)^3y_j^3}
+\end{equation}
+and
+\begin{equation}
+  b_i^{(n)}:=\mathbf v_i+2e\rho_{n-1}\mathbb U_i(n-1)^2
+\end{equation}
+in terms of which
+\begin{equation}
+  \mathbb U=A^{-1}b^{(n)}
+  .
+\end{equation}
+Finally, we compute $\rho_n$ using the second of\-~(\ref{iteration_uk}):
+\begin{equation}
+  \rho_n=\frac{2e}{\hat S_n(0)}
+  ,\quad
+  S_n(0)\approx\hat v(0)
+  -\frac1{2(2\pi)^3}\sum_{j=1}^Nw_j\frac{(1-y_j)\hat u_n(\frac{1-y_j}{y_j})\Eta(0,\frac{1-y_j}{y_j})}{y_j^3}
+  .
+\end{equation}
+
+\section{Potentials}\label{sec:potentials}
+\indent
+In this section, we describe the potentials available in {\tt simplesolv}, and provide documentation to add custom potentials to {\tt simplesolv}.
+\bigskip
+
+\subsection{Built-in potentials}
+\subsubsection{{\tt exp}}\label{sec:exp}
+\indent
+In $x$-space,
+\begin{equation}
+  v(|x|)=ae^{-|x|}
+  .
+\end{equation}
+The constant $a$ is specified through the {\tt v\_a} parameter, and can be any real number.
+Note that $v\geqslant 0$ if and only if $a\geqslant 0$.
+This is the potential that is selected by default.
+\bigskip
+
+\point
+In Fourier space,
+\begin{equation}
+  \hat v(|k|)=\int dx\ e^{ikx}v(|x|)
+  =\frac{8\pi a}{(1+k^2)^2}
+  .
+\end{equation}
+In particular, $v$ is of positive type (that is, $\hat v\geqslant 0$) if and only if $a\geqslant 0$.
+\bigskip
+
+\point
+The zero energy scattering solution, that is, the solution of
+\begin{equation}
+  (-\Delta+v)\psi=0
+  ,\quad
+  \lim_{|x|\to\infty}\psi=1
+\end{equation}
+is
+\begin{equation}
+  \psi(r)=\frac{cI_0(2\sqrt a e^{-\frac r2})+2K_0(2\sqrt a e^{-\frac r2})}r
+  ,\quad
+  c=-\frac{2K_0(2\sqrt a)}{I_0(2\sqrt a)}
+\end{equation}
+where $I_0$ and $K_0$ are the modified Bessel functions, and the square root is taken with a branch cut on the negative imaginary axis.
+In other words, if $a<0$, $\sqrt a=i\sqrt{|a|}$ and\-~\dlmfcite{(10.27.6),(10.27.9),(10.27.10)}{1.1.3}
+\begin{equation}
+  I_0(ix)=J_0(x)
+  ,\quad
+  K_0(ix)=-\frac\pi2(Y_0(x)+iJ_0(x))
+\end{equation}
+where $J_0$ and $Y_0$ are the Bessel functions, so
+\begin{equation}
+  \psi(r)=\pi\frac{c'J_0(2\sqrt{|a|} e^{-\frac r2})-Y_0(2\sqrt{|a|} e^{-\frac r2})}r
+  ,\quad
+  c'=\frac{Y_0(2\sqrt{|a|})}{J_0(2\sqrt{|a|})}
+  .
+\end{equation}
+The scattering length is\-~\dlmfcite{(10.25.2),(10.31.2)}{1.1.3}
+\begin{equation}
+  \mathfrak a_0=\lim_{r\to\infty}r(1-\psi(r))
+  =\log(a)+2\gamma+\frac{2K_0(2\sqrt a)}{I_0(2\sqrt a)}
+\end{equation}
+where $\gamma$ is the Euler constant, which, for $a<0$, is
+\begin{equation}
+  \mathfrak a_0=\log|a|+2\gamma-\frac{\pi Y_0(2\sqrt{|a|})}{J_0(2\sqrt{|a|})}
+  .
+\end{equation}
+
+\subsubsection{{\tt tent}}\label{sec:tent}
+\indent
+In $x$-space,
+\begin{equation}
+  v(|x|)=\mathds 1_{|x|<b}a\frac{2\pi}3\left(1-\frac{|x|}b\right)^2\left(\frac{|x|}b+2\right)
+  .
+\end{equation}
+The constants $a$ and $b$ are specified through the {\tt v\_a} and {\tt v\_b} parameters, and $a\in\mathbb R$, $b>0$.
+Note that $v\geqslant 0$ if and only if $a\geqslant 0$.
+Note that
+\begin{equation}
+  v(|x|)=\frac a{b^3}\mathds 1_{|x|<\frac b2}\ast\mathds 1_{|x|<\frac b2}
+\end{equation}
+(this is more easily checked in Fourier space than explicitly computing the convolution).
+\bigskip
+
+\indent
+In Fourier space,
+\begin{equation}
+  \hat v(|k|)=\int dx\ e^{ikx}v(|x|)=a\frac{b^3}8\left(4\pi\frac{\sin(\frac{|k|b}2)-\frac{|k|b}2\cos(\frac{|k|b}2)}{\frac{(|k|b)^3}8}\right)^2
+  .
+\end{equation}
+Note that this is $\frac a{b^3}(\hat\mathds 1_{|x|<\frac b2})^2$.
+In particular, $v$ is of positive type (that is, $\hat v\geqslant 0$) if and only if $a\geqslant 0$.
+
+\subsubsection{{\tt expcry}}\label{sec:expcry}
+\indent
+In $x$-space
+\begin{equation}
+  v(|x|)=e^{-|x|}-ae^{-b|x|}
+  .
+\end{equation}
+The constants $a$ and $b$ are specified through the {\tt v\_a} and {\tt v\_b} parameters, and $a\in\mathbb R$, $b>0$.
+Note that $v\geqslant 0$ if and only if $a\leqslant 1$ and $b\geqslant 1$.
+\bigskip
+
+\indent
+In Fourier space
+\begin{equation}
+  \hat v(|k|)=\int dx\ e^{ikx}v(|x|)
+  =8\pi\left(\frac1{(1+k^2)^2}-\frac{ab}{(b^2+k^2)^2}\right)
+\end{equation}
+In particular, $\hat v$ is of positive type (that is, $\hat v\geqslant 0$) if and only if $ab\leqslant 1$, $a\leqslant b$ and $a\leqslant b^3$.
+If $a\leqslant 1$, $b\geqslant 1$ and $ab>1$, then $\hat v$ has a unique minimum at
+\begin{equation}
+  |k_*|=\sqrt{\frac{b^2-(ab)^{\frac13}}{(ab)^{\frac13}-1}}
+  ,\quad
+  \hat v(|k_*|)=
+  -8\pi\frac{((ab)^{\frac13}-1)^3}{b^2-1}
+  .
+\end{equation}
+
+\subsubsection{{\tt npt}}\label{sec:npt}
+\indent
+In $x$-space
+\begin{equation}
+  v(|x|)=x^2e^{-|x|}
+  .
+\end{equation}
+Note that $v\geqslant 0$.
+\bigskip
+
+\indent
+In Fourier space
+\begin{equation}
+  \hat v(|k|)=\int dx\ e^{ikx}v(|x|)
+  =96\pi\frac{1-k^2}{(1+k^2)^4}
+  .
+\end{equation}
+In particular, $v$ is not of positive type (that is, $\hat v$ is not $\geqslant 0$).
+
+\subsubsection{{\tt alg}}\label{sec:alg}
+\indent
+In $x$-space
+\begin{equation}
+  v(|x|)=\frac1{1+\frac14|x|^4}
+  .
+\end{equation}
+Note that $v\geqslant 0$.
+\bigskip
+
+\indent
+In Fourier space
+\begin{equation}
+  \hat v(|k|)=\int dx\ e^{ikx}v(|x|)
+  =4\pi^2 e^{-|k|}\frac{\sin(|k|)}{|k|}
+  .
+\end{equation}
+In particular, $v$ is not of positive type (that is, $\hat v$ is not $\geqslant 0$).
+
+\subsubsection{{\tt algwell}}\label{sec:algwell}
+\indent
+In $x$-space
+\begin{equation}
+  v(|x|)=\frac{1+a|x|^4}{(1+|x|^2)^4}
+  .
+\end{equation}
+The constant $a$ can be set using the parameter {\tt v\_a} and can be any real number.
+Note that $v\geqslant 0$ if and only if $a\geqslant 0$.
+If $a>8$, this potential has a local minimum at $|x_-|$ and a local maximum at $|x_+|$:
+\begin{equation}
+  |x_\pm|=\sqrt{\frac12(1\pm\sqrt{1-8a^{-1}})}
+  .
+\end{equation}
+\bigskip
+
+\indent
+In Fourier space
+\begin{equation}
+  \hat v(|k|)=\int dx\ e^{ikx}v(|x|)
+  =\frac{\pi^2}{24}e^{-|k|}(a(k^2-9|k|+15)+k^2+3|k|+3)
+  .
+\end{equation}
+In particular, $v$ is of positive type (that is, $\hat v$ is not $\geqslant 0$) if and only if
+\begin{equation}
+  -\frac15\leqslant a\leqslant 3+\frac87\sqrt7\approx6.02
+  .
+\end{equation}
+
+\subsubsection{{\tt exact}}\label{sec:exact}
+\indent
+In $x$-space
+\begin{equation}
+  v(|x|)=
+  \frac{12c( x^6b^6(2e-b^2) +b^4x^4(9e-7b^2) +4b^2x^2(3e-2b^2) +(5e+16b^2))}{(1+b^2x^2)^2(4+b^2x^2)^2((1+b^2x^2)^2-c)}
+\end{equation}
+The constants $a,b,c,e$ can be set using the parameters {\tt v\_a}, {\tt v\_b}, {\tt v\_c}, {\tt v\_e}, and $a,e\in\mathbb R$, $b\neq0$, $c>0$, $c\neq9$.
+Note that $v\geqslant 0$ if and only if\-~\cite{CJL20b}
+\begin{equation}
+  b>0
+  ,\quad
+  0<c<1
+  ,\quad
+  \frac e{b^2}\geqslant\frac{-263+23\sqrt{161}}{48}\approx 0.601
+  .
+\end{equation}
+With this potential, the Simple equation has an exact solution:
+\begin{equation}
+  u=\frac c{(1+b^2x^2)^2}
+  ,\quad
+  \rho=\frac{b^3}{c\pi^2}
+  .
+\end{equation}
+\bigskip
+
+\indent
+In Fourier space $\hat v(|k|)=\int dx\ e^{ikx}v(x)$ is
+\begin{equation}
+  \begin{largearray}
+    \hat v(|k|)=\frac{48 c\pi^2}{|k|}
+    \left(
+      \frac{18+3\sqrt c-(4-3\alpha)c-(1-2\alpha)c^{\frac32}}{4(3+\sqrt c)^2c^{\frac32}}
+      e^{-\sqrt{1-\sqrt c}\frac{|k|}b}
+      \right.\\[0.5cm]\hfill\left.
+      +
+      \frac{-18+3\sqrt c+(4-3\alpha)c-(1-2\alpha)c^{\frac32}}{4(3-\sqrt c)^2c^{\frac32}}
+      e^{-\sqrt{1+\sqrt c}\frac{|k|}b}
+      +\frac{1-\frac{|k|}b}{2c}
+      e^{-\frac{|k|}b}
+      -\frac{\alpha(3(9-c)\frac{|k|}b+8c)}{8(9-c)^2}
+      e^{-2\frac{|k|}b}
+    \right)
+  \end{largearray}
+\end{equation}
+with $\alpha:=\frac e{b^2}$.
+
+\subsection{Programming custom potentials}
+\indent
+In this section, we provide documentation for programming custom potentials.
+\bigskip
+
+\indent
+The potentials are implemented in the file `{\tt \$SIMPLESOLV/src/potentials.jl}', and consist of two functions, one specifying the potential in Fourier space (the ``potential function''), and the other returning an approximate value for the scattering length (the ``scatterin glength function'') (as is explained below, a precise value of the scattering length is not actually needed).
+For instance, the potential \refname{sec:exp}{{\tt exp}} has two functions: {\tt v\_exp} and {\tt a0\_exp}.
+The potential function should take the following arguments:
+\begin{itemize}
+  \item {\tt k} ({\tt Float64}): the Fourier momentum
+  \item and any parameters that the potential depends on, such as $a$ in \refname{sec:exp}{{\tt exp}} (can be of any type, provided the appropriate changes are made to {\tt main.jl} as explained below)
+\end{itemize}
+and it must return a {\tt Float64}: the value of $\hat v$ at {\tt k}.
+The scattering length function takes the same parameters as an input, and returns a {\tt Float64}: the approximate value for the scattering length.
+\bigskip
+
+\indent
+In addition, the potential must be linked in `{\tt \$SIMPLESOLV/src/main.jl}'.
+In that file, the potential is read from the command line option {\tt U}.
+The relevant code is in lines 197-222.
+To add a new potential, add
+\begin{code}
+  elseif potential=="{\rm \{name of potential\}}"\par
+  \indent  v=k->{\rm \{potential function\}}(k,v\_param\_a,v\_param\_b,{\rm...})\par
+  \indent  a0={\rm \{scattering length function\}}(v\_param\_a,v\_param\_b,{\rm...})\par
+\end{code}
+where the number of {\tt v\_param} entries should be the number of parameters of the potential.
+The parameters that are currently read from the parameters list are {\tt a}, {\tt b}, {\tt c} and {\tt e}.
+To add a parameter, it must first be declared and initialized after line 35, and code to read it should be added after line 172:
+\begin{code}
+  elseif lhs="v\_{\rm\{name of parameter\}}"\par
+  \indent  v\_param\_{\rm\{name of parameter\}}=parse(Float64,rhs)
+\end{code}
+If the new parameter has a type other than {\tt Float64}, this should be changed in the {\tt parse} function, and in the initialization.
+\bigskip
+
+\indent
+The approximation of the scattering length is only used to initialize the Newton algorithm for \refname{sec:easyeq}{{\tt easyeq}}, so it is not important that it be exact.
+In fact, some of the built-in potentials set the scattering length to 1, when it has proved too difficult to compute it exactly.
+
+\appendix
+
+\section{Chebyshev polynomial expansion}\label{appChebyshev}
+\indent
+In this appendix, we compute the error of Chebyshev polynomial expansions of regular functions.
+Specifically, we will consider class-$s$ Gevrey functions (which is a generalization of the notion of analyticity: analytic functions are Gevrey functions with $s=1$).
+A class-$s$ Gevrey function on $[-1,1]$ is a $\mathcal C^\infty$ function that satisfies, $\forall n\in\mathbb N$,
+\begin{equation}
+  \mathop{\mathrm{sup}}_{x\in[-1,1]}\left|\frac{d^nf(x)}{dx^n}\right|\leqslant C_0C^n(n!)^s.
+\end{equation}
+Formally, the Chebyshev polynomial expansion of $f$ is
+\begin{equation}
+  f(x)=\frac{c_0}2+\sum_{j=1}^\infty c_j T_j(x)
+  \label{eqcheby}
+\end{equation}
+where $T_j$ is the $j$-th Chebyshev polynomial:
+\begin{equation}
+  T_j(x):=\cos(j\arccos(x))
+\end{equation}
+and
+\begin{equation}
+  c_j:=\frac2\pi\int_0^\pi d\theta\ f(\cos\theta)\cos(j\theta)
+  .
+\end{equation}
+
+\bigskip
+
+\theo{Lemma}\label{lemmaChebyshev}
+  Let $f$ be a class-$s$ Gevrey function on $[-1,1]$ with $s\in\mathbb N\setminus\{0\}$.
+  There exist $b_0,b>0$, which are independent of $s$, such that the coefficients $c_j$ of the Chebyshev polynomial expansion are bounded as
+  \begin{equation}
+    c_j\leqslant b_0e^{-bj^{\frac1s}}
+    .
+    \label{boundchebyshevcoef}
+  \end{equation}
+  In particular, the Chebyshev polynomial expansion is absolutely convergent pointwise (and, therefore in any $L_p$ norm), and, for every $N\geqslant 1$,
+  \begin{equation}
+    \left|f(x)-\frac{c_0}2-\sum_{j=1}^N c_j T_j(x)\right|
+    \leqslant
+    \frac{b_0}{1-e^{-b}}(s-1)!(N^{\frac1s}+b^{-1})^{s-1}e^{-bN^{\frac1s}}
+    .
+    \label{boundchebyshev}
+  \end{equation}
+\endtheo
+\bigskip
+
+\indent\underline{Proof}:
+  Note that~\-(\ref{eqcheby}) is nothing other than the Fourier cosine series expansion of $F(\theta):=f(\cos(\theta))$, which is an even, periodic, class-$s$ Gevrey function on $[-\pi,\pi]$, whose $j$-th Fourier coefficient for $j\in\mathbb Z$ is equal to $\frac12c_{|j|}$.
+  The bound\-~(\ref{boundchebyshevcoef}) follows from a well-known estimate of the decay of Fourier coefficients of class-$s$ Gevrey functions (see e.g.~\-\cite[Theorem~\-3.3]{Ta87}).
+  The bound\-~(\ref{boundchebyshev}) then follows from $|T_j(x)|\leqslant 1$ and lemma\-~\ref{lemma:ineqgevrey} below.
+\qed
+\bigskip
+
+\theo{Lemma}\label{lemma:ineqgevrey}
+  Given $b>0$ and two integers $N,s>0$,
+  \nopagebreakaftereq
+  \begin{equation}
+    \sum_{j=N}^\infty e^{-bj^{\frac1s}}
+    \leqslant
+    (s-1)!\left(N^{\frac1s}+\frac1{1-e^{-b}}\right)^{s-1}\frac{e^{-bN^{\frac1s}}}{1-e^{-b}}
+  \end{equation}
+\endtheo
+\restorepagebreakaftereq
+\bigskip
+
+\indent\underline{Proof}:
+  If $\nu_{N,s}^s:=\lfloor N^{\frac1s}\rfloor^s$ denotes the largest integer that is $\leqslant N$ and has an integer $s$-th root, then
+  \begin{equation}
+    \sum_{j=N}^\infty e^{-bj^{\frac1s}}\leqslant
+    \sum_{j=\nu_{N,s}^s}^\infty e^{-bj^{\frac1s}}\leqslant
+    \sum_{k=\nu_{N,s}}^\infty(k^s-(k-1)^s)e^{-bk}\leqslant
+    s\sum_{k=\nu_{N,s}}^\infty k^{s-1}e^{-bk}
+    .
+  \end{equation}
+  We then estimate
+  \begin{equation}
+    \sum_{k=\nu_{N,s}}^\infty k^{s-1}e^{-bk}=\frac{d^{s-1}}{d(-b)^{s-1}}\sum_{k=\nu_{N,s}}^\infty e^{-bk}
+    \leqslant (s-1)!\left(\nu_{N,s}+\frac1{1-e^{-b}}\right)^{s-1}\frac{e^{-b\nu_{N,s}}}{1-e^{-b}}
+  \end{equation}
+  which concludes the proof of the lemma.
+\qed
+
+\section{Gauss quadratures}\label{appGL}
+\indent
+Gauss quadratures are approximation schemes to compute integrals of the form
+\begin{equation}
+  \int_a^bdt\ \omega(t)f(t)
+\end{equation}
+where $\omega(x)\geqslant 0$ is one of several functions that Gauss quadratures can treat.
+The possible choices of $\omega$ and $(a,b)$ are
+\begin{itemize}
+  \item $(a,b)=(-1,1)$, $\omega(t)=1$: Gauss-Legendre quadratures.
+  \item $(a,b)=(-1,1)$, $\omega(t)=(1-t)^\alpha(1+t)^\beta$, $\alpha,\beta>-1$: Gauss-Jacobi quadratures.
+  \item $(a,b)=(-1,1)$, $\omega(t)=(1-t^2)^{-\frac12}$: Gauss-Chebyshev quadratures.
+  \item $(a,b)=(0,\infty)$, $\omega(t)=e^{-t}$: Gauss-Laguerre quadratures.
+  \item $(a,b)=(0,\infty)$, $\omega(t)=t^\alpha e^{-t}$, $\alpha>-1$: generalized Gauss-Laguerre quadratures.
+  \item $(a,b)=(0,\infty)$, $\omega(t)=e^{-t^2}$: Gauss-Hermite quadratures.
+\end{itemize}
+It is not our goal here to discuss Gauss quadratures in detail, or their relation to orthogonal polynomials.
+Instead, we will compute the error made when approximating such an integral by a Gauss quadrature.
+\bigskip
+
+\indent
+For each Gauss quadrature, the integral is approximated in the form
+\begin{equation}
+  \int_a^b dt\ \omega(t) f(t)
+  \approx
+  \sum_{i=1}^Nw_if(r_i)
+\end{equation}
+where $w_i$ are called the {\it weights}, $r_i$ are the {\it abcissa}, and $N$ is the {\it order}.
+The weights and abcissa depend on both $\omega$ and the order $N$.
+The crucial property of Gauss quadratures is that they are {\it exact} when $f$ is a polynomial of order $\leqslant 2N-1$.
+\bigskip
+
+\indent
+In this appendix, we compute the error of Gauss quadratures when used to integrate regular functions.
+Specifically, we will consider class-$s$ Gevrey functions (which is a generalization of the notion of analyticity: analytic functions are Gevrey functions with $s=1$).
+A class-$s$ Gevrey function on is a $\mathcal C^\infty$ function that satisfies, $\forall n\in\mathbb N$,
+\begin{equation}
+  \mathop{\mathrm{sup}}_{x\in[-1,1]}\left|\frac{d^nf(x)}{dx^n}\right|\leqslant C_0C^n(n!)^s.
+\end{equation}
+
+\bigskip
+
+\theo{Lemma}\label{lemmaGL}
+  Let $f$ be a class-$s$ Gevrey function with $s\in\mathbb N\setminus\{0\}$.
+  There exist $b_0,b>0$, which are independent of $s$, and $N_0>0$, which is independent of $s$ and $f$, such that, if $N\geqslant N_0$, then, denoting the Gauss weights and abcissa by $w_i$ and $r_i$,
+  \begin{equation}
+    \left|\int_a^bdt\ \omega(t)f(t)-\sum_{i=1}^Nw_if(r_i)\right|
+    \leqslant
+    b_0(s-1)!((2N-1)^{\frac1s}+b^{-1})^{s-1}e^{-b(2N-1)^{\frac1s}}\left(\int_a^b\omega(t)+\sum_{i=1}^Nw_i\right)
+    .
+  \end{equation}
+  In particular, if $f$ is analytic (i.e. $s=1$), then
+  \begin{equation}
+    \left|\int_a^bdt\ \omega(t)f(t)-\sum_{i=1}^Nw_if(r_i)\right|
+    \leqslant b_0e^{-b(2N-1)}\left(\int_a^b\omega(t)+\sum_{i=1}^Nw_i\right)
+    .
+  \end{equation}
+\endtheo
+\bigskip
+
+\indent\underline{Proof}:
+  We expand $f$ into Chebyshev polynomials as in\-~(\ref{eqcheby}):
+  \begin{equation}
+    f(x)=\frac{c_0}2+\sum_{j=1}^{\infty}c_jT_j(x)
+  \end{equation}
+  Let
+  \begin{equation}
+    g(x):=f(x)-\frac{c_0}2-\sum_{j=1}^{2N-1} c_j T_j(x)
+    .
+  \end{equation}
+  Since order-$N$ Gauss quadratures are exact on polynomials of order $\leqslant 2N-1$, we have
+  \begin{equation}
+    \int_a^bdt\ \omega(t)f(t)-\sum_{i=1}^Nw_if(r_i)
+    =
+    \int_a^bdt\ \omega(t)g(t)-\sum_{i=1}^Nw_ig(r_i)
+  \end{equation}
+  and, by lemma\-~\ref{lemmaChebyshev},
+  \begin{equation}
+    |g(x)|\leqslant
+    (\mathrm{const}.)(s-1)!((2N-1)^{\frac1s}+b^{-1})^{s-1}e^{-b(2N-1)^{\frac1s}}
+    .
+    \label{boundchebyshev}
+  \end{equation}
+\qed
+
+\section{Bipolar coordinates}\label{appendix:bipolar}
+\indent
+Bipolar coordinates are very useful for computing convolutions of radial functions in three dimensions.
+\bigskip
+
+\theo{Lemma}\label{lemma:bipolar}
+  For $y\in\mathbb R^3$,
+  \nopagebreakaftereq
+  \begin{equation}
+    \int_{\mathbb R^3}dx\ f(|x|,|x-y|)
+    =
+    \frac{2\pi}{|y|}\int_0^\infty dt\int_{||y|-t|}^{|y|+t}ds\ stf(s,t)
+  \end{equation}
+\endtheo
+\restorepagebreakaftereq
+
+\indent\underline{Proof}:
+  Without loss of generality, we assume that $y=(0,0,a)$ with $a>0$. We first change to cylindrical coordinates: $(\rho,\theta,x_3)$:
+  \begin{equation}
+    \int_{\mathbb R^3}dx\ f(|x|,|x-y|)
+    =
+    2\pi\int_0^\infty d\rho\int_{-\infty}^\infty dx_3\ \rho f(|(\rho,0,x_3)|,|(\rho,0,x_3-a)|)
+    .
+  \end{equation}
+  Next, we change variables to
+  \begin{equation}
+    s=|(\rho,0,x_3)|
+    ,\quad
+    t=|(\rho,0,x_3-a)|
+    .
+  \end{equation}
+  The inverse of this transformation is
+  \begin{equation}
+    x_3=\pm\frac{a^2+s^2-t^2}{2a}
+    ,\quad
+    \rho=\frac1{2a}\sqrt{4a^2s^2-(a^2+s^2-t^2)^2}
+    \label{inverse}
+  \end{equation}
+  and its Jacobian is
+  \begin{equation}
+    \frac{2ts}{\sqrt{-2a^4+2a^2(t^2+s^2)-(t^2-s^2)^2}}
+    =\frac{ts}{\rho a}
+    .
+    \label{jacobian}
+  \end{equation}
+\qed
+\bigskip
+
+\indent
+The following is a generalization of the previous lemma to functions of four variables.
+\bigskip
+
+\theo{Lemma}\label{lemma:bipolar2}
+  For $y\in\mathbb R^3$,
+  \begin{equation}
+    \begin{largearray}
+      \int_{\mathbb R^6}dxdx'\ f(|x|,|x-y|,|x'|,|x'-y|,|x-x'|)
+      \\[0.5cm]\hfill
+      =
+      \frac{2\pi}{|y|^2}\int_0^\infty dt\int_{||y|-t|}^{|y|+t}ds\int_0^\infty dt'\int_{||y|-t'|}^{|y|+t'}ds'\int_0^{2\pi}d\theta\ sts't'f(s,t,s',t',\xi(s,t,s',t',\theta,|y|))
+    \end{largearray}
+  \end{equation}
+  with
+  \nopagebreakaftereq
+  \begin{equation}
+    \begin{largearray}
+      \xi(s,t,s',t',\theta,|y|)
+      :=
+      \frac1{\sqrt 2|y|}
+      \Big(
+	|y|^2(s^2+t^2+(s')^2+(t')^2)-|y|^4-(s^2-t^2)((s')^2-(t')^2)
+	\\[0.5cm]\hfill
+	-
+	\sqrt{(4|y|^2s^2-(|y|^2+s^2-t^2)^2)(4|y|^2(s')^2-(|y|^2+(s')^2-(t')^2)^2)}\cos\theta
+      \Big)^{\frac12}
+      .
+    \end{largearray}
+  \end{equation}
+\endtheo
+\restorepagebreakaftereq
+
+\indent\underline{Proof}:
+  Without loss of generality, we assume that $y=(0,0,a)$ with $a>0$. We first change to cylindrical coordinates: $(\rho,\theta,x_3;\rho',\theta',x_3')$:
+  \begin{equation}
+    \begin{largearray}
+      \int_{\mathbb R^6}dxdx'\ f(|x|,|x-y|,|x'|,|x'-y|,|x-x'|)
+      \\[0.5cm]\hfill
+      =
+      2\pi\int_0^\infty d\rho\int_{-\infty}^\infty dx_3
+      \int_0^\infty d\rho'\int_{-\infty}^\infty dx_3'\int_0^{2\pi}d\theta'
+      \ \rho\rho'
+      f(
+	s,t,s',t',
+	|(\rho-\rho'\cos\theta',\rho'\sin\theta',x_3-x_3')|
+      )
+      .
+    \end{largearray}
+  \end{equation}
+  where
+  \begin{equation}
+    s:=|(\rho,0,x_3)|
+    ,\quad
+    t:=|(\rho,0,x_3-a)|
+    ,\quad
+    s':=|(\rho'\cos\theta',\rho'\sin\theta',x_3')|,
+    t':=|(\rho'\cos\theta',\rho'\sin\theta',x_3'-a)|
+  \end{equation}
+  Next, we change variables to $(s,t,s',t',\theta')$.
+  The Jacobian of this transformation is, by\-~(\ref{jacobian}),
+  \begin{equation}
+    \frac{tst's'}{\rho\rho' a^2}
+    .
+  \end{equation}
+  Furthermore, by\-~(\ref{inverse}),
+  \begin{equation}
+    x_3-x_3'=\frac{s^2-t^2-(s')^2+(t')^2}{2a}
+  \end{equation}
+  and
+  \begin{equation}
+    \rho=\frac{\sqrt{4a^2s^2-(a^2+s^2-t^2)^2}}{2a}
+    ,\quad
+    \rho'=\frac{\sqrt{4a^2(s')^2-(a^2+(s')^2-(t')^2)^2}}{2a}
+  \end{equation}
+  so
+  \begin{equation}
+    \begin{largearray}
+      |(\rho-\rho'\cos\theta',\rho'\sin\theta',x_3-x_3')|
+      =
+      \frac1{2a}
+      \left(
+	4a^2s^2-(a^2+s^2-t^2)^2
+	+
+	4a^2(s')^2-(a^2+(s')^2-(t')^2)^2
+	\right.\\[0.5cm]\hfill\left.
+	-
+	2\sqrt{(4a^2s^2-(a^2+s^2-t^2)^2)(4a^2(s')^2-(a^2+(s')^2-(t')^2)^2)}\cos\theta'
+	+
+	(s^2-t^2-(s')^2+(t')^2)^2
+      \right)^{\frac12}
+    \end{largearray}
+  \end{equation}
+  and
+  \begin{equation}
+    \begin{largearray}
+      |(\rho-\rho'\cos\theta',\rho'\sin\theta',x_3-x_3')|
+      =
+      \frac1{\sqrt 2a}
+      \Big(
+	a^2(s^2+t^2+(s')^2+(t')^2)-a^4-(s^2-t^2)((s')^2-(t')^2)
+	\\[0.5cm]\hfill
+	-
+	\sqrt{(4a^2s^2-(a^2+s^2-t^2)^2)(4a^2(s')^2-(a^2+(s')^2-(t')^2)^2)}\cos\theta'
+      \Big)^{\frac12}
+      .
+    \end{largearray}
+  \end{equation}
+\qed
+
+\section{Hann windowing}\label{appendix:hann}
+\indent
+In this appendix, we discuss the use of Hann windows to compute Fourier transforms.
+Consider the Fourier transform
+\begin{equation}
+  \hat f(k)=\int dx\ e^{ikx}f(x)
+  .
+\end{equation}
+Evaluating this integral numerically can be tricky, especially at high $|k|$, because of the rapid oscillations at large $|x|$.
+A trick to palliate such a problem is to multiply $f$ by a {\it window function} $h_L$, which cuts off distances of order $L$.
+We then compute, instead of $\hat f$,
+\begin{equation}
+  \tilde f(k)=\int dx\ e^{ikx}h_L(x)f(x)
+  .
+\end{equation}
+We can then evaluate $\tilde f$ using standard numerical techniques, such as Gauss quadratures (see appendix\-~\ref{appGL}), without issues at large $|x|$.
+However, in doing so, we will make an error in the Fourier transform.
+To quantify this error, note that
+\begin{equation}
+  \tilde f(k)
+  =\hat h_L\hat\ast\hat f(k)
+  \equiv\int\frac{dq}{(2\pi)^d}\ \hat h_L(q)\hat f(k-q)
+\end{equation}
+so if we choose $h_L$ in such a way that $\hat h_L$ is peaked around the origin, then $\tilde f$ will not differ too much from $\hat f$:
+\begin{equation}
+  \hat f(k)-\tilde f(k)
+  =((2\pi)^d\delta(k)-\hat h_L)\hat\ast\hat f(k)
+  .
+\end{equation}
+\bigskip
+
+\indent
+The {\it Hann} window is defined as
+\begin{equation}
+  h_L(x)=\cos^2({\textstyle\frac{\pi|x|}L})\mathds 1_{|x|<\frac L2}
+\end{equation}
+whose Fourier transform is, in $d=3$,
+\begin{equation}
+  \hat h_L(k)=
+  \frac{4\pi^3 L^2}{|k|}\frac{((|k|L)^3-4|k|L\pi^2)\cos(\frac{|k|L}2)-2(3(|k|L)^2-4\pi^2)\sin(\frac{|k|L}2)}{((|k|L)^3-4|k|L\pi^2)^2}
+\end{equation}
+which decays at large $|k|L$ as
+\begin{equation}
+  \hat h_L(k)\sim\frac{4\pi^3}{|k|^4L}\cos({\textstyle\frac{|k|L}2})
+  .
+\end{equation}
+
+
+
+\vfill
+\eject
+
+\begin{thebibliography}{WWW99}
+\small
+\IfFileExists{bibliography/bibliography.tex}{\input bibliography/bibliography.tex}{}
+\end{thebibliography}
+
+
+\end{document}