I'm looking for closed form of points where $\operatorname{Erfc}(-x)$ intersects with $a \cdot e^{-x^2}$ for $a$ being some fixed non-zero constant.
I guess I can use the exponential approximation for $\operatorname{Erfc}$. Then I probably can find closed form with help of Lambert W function. But I wonder if solution can be found without any use of approximations for $\operatorname{Erfc}$.
Say that you want to find the inverse of $$a=e^{x^2} (1+\text{erf}(x))$$ which, I suppose, does not show any closed form solution.
Better is to consider instead the problem of the inverse of $$b=\log\Big(e^{x^2} (1+\text{erf}(x))\Big) \quad \text{with} \quad b=\log(a)$$
What we can notice is that $\large\color{red}{x_0=\sqrt b}$ is a decent approximation. For example, for $b=4$, the "exact" solution is $1.81986$. Moreover $$0\leq\log \left(e^b \left(1+\text{erf}\left(\sqrt{b}\right)\right)\right)-b < \log(2)$$
By Darboux theorem, we also know that, $x_0$ being an over estimate, the solution by Newton method will converge without any overshoot. For the given example, $x_1=1.82774$, $x_2=1.81988$.
This means that we can have $\color{red}{\text{explicit}}$ approximate solutions using the first iterate of Newton-like methods of order $n$ starting at $x_0=\sqrt b$.
A few numbers for illustration $$\left( \begin{array}{cccccc} b & x_1^{(2)} & x_1^{(3)}& x_1^{(4)}& x_1^{(5)} & \text{solution}\\ 0.5 & 0.42132 & 0.38998 & 0.38861 & 0.38872 & 0.38872 \\ 1.0 & 0.72532 & 0.69733 & 0.69579 & 0.69586 & 0.69586 \\ 1.5 & 0.97265 & 0.95039 & 0.94912 & 0.94914 & 0.94915 \\ 2.0 & 1.18365 & 1.16622 & 1.16525 & 1.16525 & 1.16525 \\ 2.5 & 1.36912 & 1.35533 & 1.35461 & 1.35460 & 1.35460\\ 3.0 & 1.53563 & 1.52454 & 1.52400 & 1.52399 & 1.52399 \\ 3.5 & 1.68751 & 1.67842 & 1.67802 & 1.67801 & 1.67801 \\ 4.0 & 1.82774 & 1.82018 & 1.81987 & 1.81986 & 1.81986 \\ 4.5 & 1.95850 & 1.95211 & 1.95186 & 1.95186 & 1.95186 \\ 5.0 & 2.08138 & 2.07589 & 2.07570 & 2.07570 & 2.07570 \\ \end{array} \right)$$
Edit
It is possible to find a better $x_0$ since the solution is such that $$\sqrt{b-\log (2)} < x < \sqrt{b}$$
So, working an initial estimate of the form $$x_0=\sqrt{b-\log (2) \,f(b)} $$ such that $f(0)=0$ and $f(\infty)=1$
Some empirical analysis leads to propose the simplistic $$f(b)=\text{erf}\left(\frac{253}{321} x^{\frac 7{10}}\right)$$ This was obtained by a quick and dirty nonlinear regression for $0\leq b \leq 5$ which correspond to an adjusted $R^2=0.99995$.
Repeating the calculations only for Newton method
$$\left( \begin{array}{cccc} b & x_0 & x_1^{(2)} & \text{solution}\\ 0.5 & 0.38512 & 0.38873 & 0.38872 \\ 1.0 & 0.70039 & 0.69587 & 0.69586 \\ 1.5 & 0.95028 & 0.94915 & 0.94915 \\ 2.0 & 1.16426 & 1.16525 & 1.16525 \\ 2.5 & 1.35300 & 1.35460 & 1.35460 \\ 3.0 & 1.52252 & 1.52399 & 1.52399 \\ 3.5 & 1.67689 & 1.67801 & 1.67801 \\ 4.0 & 1.81910 & 1.81986 & 1.81986 \\ 4.5 & 1.95136 & 1.95186 & 1.95186 \\ 5.0 & 2.07539 & 2.07570 & 2.07570 \\ \end{array} \right)$$
This answer is probably not in the spirit of what the asker requested ("without any use of approximations for $\mathrm{Erfc}$"), but it may help the next person with a similar problem and no such restriction.
A closed-form solution does not exist. However, many computational environments provide an exponentially scaled complementary error function $e^{x^2}\mathrm{erfc}(x)$, often named erfcx, that one can use to compute $a$ from $x$. Where it is not provided, one could derive it from this sample implementation in C that I posted on Stackoverflow some time ago.
Here we need to solve $\mathrm{erfcx}(-x) = a$, that is, compute $x = -\mathrm{erfcx}^{-1}(a)$. While $\mathrm{erfcx}^{-1}$ (named erfcxinv) is not typically provided as a standard function, one can easily compute it by Newton-Raphson iteration. One observes that for $x \to 0$, $\mathrm{erfcx}^{-1}(x) \approx \frac{1}{x\sqrt{\pi}}$, and that for $x \to \infty$, $\mathrm{erfcx}^{-1}(x) \approx -\sqrt{\ln\left(\frac{x}{2}\right)}$ and can pick a starting value for the iteration process accordingly. To compute x = erfcxinv (a), one proceeds as follows:
c := 1 / sqrt (pi)
IF a < 1
x := c / a
ELSE
x := -sqrt (ln (a))
FI
DO i = 0 to n
e := erfcx (x)
r := 1 / fma (x, e, -c)
x := fma ((a - e) / 2, r, x)
OD
In the above pseudocode, fma (a,b,c) := a * b + c, with many computational environments supplying a numerically beneficial operation called fused multiply-add under that name. For IEEE-754 double-precision computation, one might choose $n = 5$ to compute results accurate to machine precision.
