  We study the fully explicit numerical approximation of a semilinear elliptic
boundary value model problem, which features a monomial reaction and analytic
forcing, in a bounded polygon $\Omega\subset\mathbb{R}^2$ with a finite number
of straight edges. In particular, we analyze the convergence of $hp$-type
iterative linearized Galerkin ($hp$-ILG) solvers. Our convergence analysis is
carried out for conforming $hp$-finite element (FE) Galerkin discretizations on
sequences of regular, simplicial partitions of $\Omega$, with geometric corner
refinement, with polynomial degrees increasing in sync with the geometric mesh
refinement towards the corners of $\Omega$. For a sequence of discrete
solutions generated by the ILG solver, with a stopping criterion that is
consistent with the exponential convergence of the exact $hp$-FE Galerkin
solution, we prove exponential convergence in $\mathrm{H}^1(\Omega)$ to the
unique weak solution of the boundary value problem. Numerical experiments
illustrate the exponential convergence of the numerical approximations obtained
from the proposed scheme in terms of the number of degrees of freedom as well
as of the computational complexity involved.
