I'm trying to convolve two functions $f$ and $g$. $$f(x) = e^{-\frac{{(x-p_2)}^2}{2 q_2^2}}$$ $$g(x) = \left(i_1 e^{-\frac{(a-x)^2}{2 \sigma ^2}}+j_1 e^{-\frac{(b-x)^2}{2 \sigma ^2}}\right) \left(i_0 e^{-\frac{\left(a+p_1-x\right){}^2}{2 \left(q_1^2+\sigma ^2\right)}}+j_0 e^{-\frac{\left(b+p_1-x\right){}^2}{2 \left(q_1^2+\sigma ^2\right)}}\right)$$
The convolution is then this: $$ \int_{-\infty }^{\infty } e^{-\frac{{(x-p_2)}^2}{2 q_2^2}} \left(\left(i_1 e^{-\frac{(a-x)^2}{2 \sigma ^2}}+j_1 e^{-\frac{(b-x)^2}{2 \sigma ^2}}\right) \left(i_0 e^{-\frac{\left(a+p_1-x\right){}^2}{2 \left(q_1^2+\sigma ^2\right)}}+j_0 e^{-\frac{\left(b+p_1-x\right){}^2}{2 \left(q_1^2+\sigma ^2\right)}}\right)\text{/.}\, x\to z-x\right) \, dx $$
These are the assumptions: $\left(\sigma |a|b|i_0|i_1|j_0|j_1|p_0|p_1|q_0|q_1\right)\in \mathbb{R}\land \sigma >0\land q_1>0\land q_2>0$
Fiddling with Mathematica and intuition I've got to this: $$\sqrt{2 \pi } q_2 \sigma \sqrt{\frac{q_1^2+\sigma ^2}{2 q_2^2 \sigma ^2+q_1^2 \left(q_2^2+\sigma ^2\right)+\sigma ^4}} \times \left(i_0 e^{-\frac{\left(x-a-p_1-p_2\right){}^2}{2 \left(\sigma ^2+q_1^2+q_2^2\right)}}+j_0 e^{-\frac{\left(x-b-p_1-p_2\right){}^2}{2 \left(\sigma ^2+q_1^2+q_2^2\right)}}\right) \left(i_1 e^{-\frac{\left(x-a-p_2\right){}^2}{2 \left(\sigma ^2+q_2^2\right)}}+j_1 e^{-\frac{\left(x-b-p_2\right){}^2}{2 \left(\sigma ^2+q_2^2\right)}}\right)$$
It's seemingly very close, however it's not correct.
It doesn't quite match the correct curve:

Two questions:
- Any idea what's wrong?
- What's the correct way to do this, without relying on Mathematica and intuition?
What's wrong is that you did the convolution wrong. Your g(x) has two factors which I'll call g1(x) and g2(x). What you appear to have calculated is $(f(x)*g1(x))\times(f(x)*g2(x))$ where * is the convolution operation. That's just not correct. What you need to do is:
Now you will have four Gaussians and you can perform the integral on each one separately. The necessary algebra is high-school level, but there will be a lot of it.