discretization  are solved by means of a direct solver; in this way, this approach differs from inexact Newton methods as discussed, e.g., in [@ErVo13]; see also [@5].

-   The parameters are chosen to be $\gamma=0.5$ and $\tau=0.1$ in Algorithm \[al:zs\], $\vartheta=0.5$ in , $\zeta=0.6$ in , and $\theta=0.5$ in . Unless stated differently, the procedure is initiated with a uniform initial mesh ${\mathcal{T}}$ consisting of 10 elements, and a polynomial degree distribution $\bm p=(1,\ldots 1)$.

\[ex:slin0\] We begin by looking at the (non-singularly perturbed) *Bratu* equation given by $$\label{eq:Bratu}
-u''=\exp(u+1)\quad\text{on }\Omega=(0,1),\qquad u(0)=u(1)=0.$$ This problem features two positive and concave solutions, which we aim to find by means of the fully adaptive Newton-Galerkin procedure introduced above (cf. also [@AmreinWihler:14 Example 3.4]). To this end, we let the initial guesses to be ${u^{\rm hp}}_0\equiv0$ and ${u^{\rm hp}}_0(x)=\pi_1(10x(1-x))$, where $\pi_1$ denotes the piecewise linear interpolation operator on the initial mesh. For these choices, Algorithm \[al:full\] generates the two solutions in bold and dashed lines, respectively, in Figure \[fig:slin0\] (left). Since both solutions are smooth, the $hp$-adaptive testing strategy  dictates pure $p$-refinements in all steps as expected, and, consequently, the residuals decay at an exponential rate with respect to the number of degrees of 