Let $X$ be the constrained random walk on ${mathbb Z}_+^d$ representing the queue lengths of a stable Jackson network and $x$ its initial position. Let $tau_n$ be the first time the sum of the components of $X$ equals $n$. $p_n doteq P_x(tau_n < tau_0)$ is a key performance measure for the queueing system represented by $X$, stability implies $p_nrightarrow 0$ exponentially. Currently the only analytic method available to approximate $p_n$ is large deviations analysis, which gives the exponential decay rate of $p_n$. Finer results are available via rare event simulation. The present article develops a new method to approximate $p_n$ and related expectations. The method has two steps: 1) with an affine transformation, move the origin onto the exit boundary of $tau_n$, take limits to remove some of the constraints on the dynamics, this yields a limit unstable constrained walk $Y$ 2) Construct a basis of harmonic functions of $Y$ and use them to apply the classical superposition principle of linear analysis. The basis functions are linear combinations of $log$-linear functions and come from solutions of harmonic systems, which are graphs whose vertices represent points on the characteristic surface of $Y$, the edges between the vertices represent conjugacy relations between the points, the loops represent membership in the boundary characteristic surfaces. Using our method we derive explicit, simple and almost exact formulas for $P_x(tau_n < tau_0)$ for $d$-tandem queues, similar to the product form formulas for the stationary distribution of $X$. The same method allows us to approximate the Balayage operator mapping $f$ to $x rightarrow {mathbb E}_x left[ f(X_{tau_n}) 1_{{tau_n < tau_0}} right]$ for a range of stable constrained random walks in $2$ dimensions. We indicate how the ideas of the paper relate to more general processes and exit boundaries.