  We present a quantum algorithm to solve systems of linear equations of the
form $A\mathbf{x}=\mathbf{b}$, where $A$ is a tridiagonal Toeplitz matrix and
$\mathbf{b}$ results from discretizing an analytic function, with a circuit
complexity of $poly(\log(\kappa), 1/\sqrt{\epsilon}, \log(N))$, where $N$
denotes the number of equations, $\epsilon$ is the accuracy, and $\kappa$ the
condition number. The \emph{repeat-until-success} algorithm has to be run
$\mathcal{O}\left(\kappa/(1-\epsilon)\right)$ times to succeed, leveraging
amplitude amplification, and sampled $\mathcal{O}(1/\epsilon^2)$ times. Thus,
the algorithm achieves an exponential improvement with respect to $N$ over
classical methods. In particular, we present efficient oracles for state
preparation, Hamiltonian simulation and a set of observables together with the
corresponding error and complexity analyses. As the main result of this work,
we show how to use Richardson extrapolation to enhance Hamiltonian simulation,
resulting in an implementation of Quantum Phase Estimation (QPE) within the
algorithm with $1/\sqrt{\epsilon}$ circuit complexity instead of $1/\epsilon$
and which can be parallelized. Furthermore, we analyze necessary conditions for
the overall algorithm to achieve an exponential speedup compared to classical
methods. Our approach is not limited to the considered setting and can be
applied to more general problems where Hamiltonian simulation is approximated
via product formulae, although, our theoretical results would need to be
extended accordingly. All the procedures presented are implemented with Qiskit
and tested for small systems using classical simulation as well as using real
quantum devices available through the IBM Quantum Experience.
