Let's propose you have the sum of two error functions:
$f(x) = \text{erf}(ax)+\text{erf(bx)}$
If you wanted to solve for x, you'd first unify them into a third function
$f(x) = \text{erf}(cx)$
Or something similar. How could you approach this problem?
$\DeclareMathOperator{erf}{erf}$ To go from your first to your second equation, we want $$ \erf(a x) + \erf(b x) = \erf(c x) \tag{1}\label{original} $$ Ideally, we could construct a function $c(a,b)$, depending only on $a$ and $b$, that satisfies \ref{original}. Let's see if that is possible. Let's relax the equation to $$ \erf(a x) + \erf(b x) = \erf(c y) $$ to which we will specialize $y = x$ later to recover solutions to \ref{original}. Let's get an idea of what we are wading into. We'll set $a = -1$, $b = -2$ and plot contours for various values of $x$ on the $cy$-plane, highlighting the points where $x = y$.
There are a few takeaways from this plot:
The last point means that a valid way forward is to solve \ref{original} for $c$ and announce that we are done. $$ c = \frac{\erf^{-1} \left( \erf(ax) + \erf(bx) \right)}{x} \text{,} $$ where we use the inverse error function, $\erf^{-1}$. If you put this into $\ref{original}$ you quickly realize that this is not responsive to your core question, which is solving $f(x) = \erf(ax) + \erf(bx)$.
Since there is no hope of an addition formula for $\erf$, and you have given no hints on the structure of $f$, we have to go to generic numerical techniques for solving. First, constrain ranges. Since the range of $\erf$ is $[-1,1]$, the right-hand side of $$ f(x) = \erf(ax) + \erf(bx) \label{eqn}\tag{2} $$ lies in $[-2,2]$, and so must the left-hand side. If you can solve $$ -2 \leq f(x) \leq 2 $$ you will get a set of intervals, each of which contains zero or more solutions in $x$. (An interval need not contain a solution since we have made no attempt to get the left-hand and right-hand sides of \ref{eqn} to agree.)
If you cannot solve that inequality, then use normal calculus optimization techniques to partition the $x$ axis into intervals on which $f(x)$ is either increasing or decreasing. For each turning point (where the graph switches from rising to falling or vice versa), check whether it is above, in, or below the band $[-2,2]$. An interval with both endpoints above, or both endpoints below, cannot contain a solution. This produces a collection of intervals in which there may be a solution to \ref{eqn}.
For each interval identified above, repeat the process for $f(x) - (\erf(ax) + \erf(bx))$ where now we keep intervals whose turning points (i.e., endpoints) have opposing sign, because the band we are interested in is $[0,0]$. Each resulting interval is guaranteed to contain an $x$ satisfying \ref{eqn}. Such an interval may contain more than one solution.
Let $f(x) = x^2$, so we study $x^2 = \erf(-x) + \erf(-2x)$. We can readily solve $-2 \leq x^2 \leq 2$ -- we immediately know that any solution satsifies $-\sqrt{2} \leq x \leq \sqrt{2}$. The first derivative of $x^2 -( \erf(-x) + \erf(-2x))$ is $2x + \frac{2 e^{-x^2}}{\sqrt{\pi }} + \frac{4 e^{-4 x^2}}{\sqrt{\pi }}$. Plotting that derivative over $[-\sqrt{2}, \sqrt{2}]$, we see
So there is one zero of that derivative. Running bisection on the interval $[-\sqrt{2}, 0]$ (on which interval the graph shows the derivative is monotonically increasing), we find the zero of the first derivative is at $x = -0.6224{\dots}$. We can check the derivative at test points on either side of this zero, or observe from the graph that $x^2 -( \erf(-x) + \erf(-2x))$ monotonically decreases to the left of this critical point and increases to the right. Thus, we get two intervals for bisection, $[-\sqrt{2}, -0.6224]$ and $[-0.6224, \sqrt{2}]$. (Pro tip: when feeding numerical endpoints to a binary search, actually use numbers very slightly past the endpoint to ensure that representation/quantization error doesn't inadvertently produce a smaller interval than intended.) For the first interval, bisection finds the solution $x = -1.3970{\dots}$; for the second, $x = 0$. (You may actually get a small deviation from zero on any particular run of a bisection. For instance, my first, coarse run gave $x = 6.6{\dots} \times 10^{-21}$. This is suspiciously close to a round number, so I re-ran it with a higher precision and accuracy goal, yielding ever more accurate approximations to zero. So either zero is a solution or it is an astoundingly good approximation of a solution.)
Let's check those. \begin{align*} (-1.3970)^2 &= 1.9516{\dots} \\ \erf(-(-1.3970)) + \erf(-2(-1.3970)) &= 1.9517{\dots} \end{align*} (Using arbitrary precision code, we can sharpen this to $x = -1.397\,045\,412\,320\,144\,210\,312\,866\,479\,771\,694$, which gives a discrepancy of $<10^{-33}$ between the left-hand and right-hand sides of \ref{eqn}.) Also,\begin{align*} (0)^2 &= 0 \\ \erf(-(0)) + \erf(-2(0)) &= 0 \end{align*} From which, we see that our second approximate solution is exact. Our method ensures that we have found all solutions, but also,
shows that we have found both solutions and that we had the correct intervals in each step.
As already said in comments and answer, you will need some numerical method for solving for $x$ the nonlinear equation $$y=\text{erf}(a x)+\text{erf}(b x)$$
I shall consider the case where $a$ and $b$ are positive (the other situations being more complex).
The first point is to notice that $y$ is bounded by $2\text{erf}(a x)$ and $2\text{erf}(b x)$.
Now, to have bounds, we could use the simple approximation $$\text{erf}(t)\approx \sqrt{1-\exp\Big(-\frac 4 {\pi}t^2 \Big)}$$ which makes the solution to be between $$\frac{\sqrt{\pi }}{2} \sqrt{-\frac{\log \left(1-\frac{y^2}{4}\right)}{a^2}}\qquad \text{and} \qquad \frac{\sqrt{\pi }}{2} \sqrt{-\frac{\log \left(1-\frac{y^2}{4}\right)}{b^2}}$$
Start at the midpoint and use Newton method. Using for example $a=0.5$, $b=1.5$ and $y=1.2345$, the iterates would be $$\left( \begin{array}{cc} n & x_n \\ 0 & 0.81835949 \\ 1 & 0.67741900 \\ 2 & 0.69292373 \\ 3 & 0.69314665 \\ 4 & 0.69314670 \end{array} \right)$$ What we could also do is to expand as Taylor series $$y=\frac{2 x (a+b)}{\sqrt{\pi }}-\frac{2 x^3 \left(a^3+b^3\right)}{3 \sqrt{\pi }}+\frac{x^5 \left(a^5+b^5\right)}{5 \sqrt{\pi }}-\frac{x^7 \left(a^7+b^7\right)}{21 \sqrt{\pi }}+O\left(x^9\right)$$ and use series reversion to get $$x=t+\alpha t^3+\beta t^5+\gamma t^7+O\left(t^9\right)\qquad \text{and} \qquad t=\frac{\sqrt{\pi } }{2 (a+b)}y$$ $$\alpha=\frac{1}{3} \left(a^2-a b+b^2\right)$$ $$\beta=\frac{1}{30} \left(7 a^4-17 a^3 b+27 a^2 b^2-17 a b^3+7 b^4\right)$$ $$\gamma=\frac{1}{630} \left(127 a^6-519 a^5 b+1191 a^4 b^2-1471 a^3 b^3+1191 a^2 b^4-519 a b^5+127 b^6\right)$$ For the worked example, this would give, as an estimate to start with, $x=0.685654$.
