Element-wise powering vs. discrete version of binomial theorem on a sum of Gaussians
([a b c] .+ [d e f]).^2 = [a+d b+e c+f].^2
                        = [(a+d)^2 (b+e)^2 (c+f)^2]
                        = [a^2+2ad+d^2 b^2+2be+e^2 c^2+2cf+f^2]

I suspect that:

The problem may be in your assumptions: &quot;We should see three peaks in this case&quot; sounds like a hypothesis that may be false.
You may have a bug in your code, which might be why &quot;We get only two peaks as a result&quot;.

Here is some code to illustrate that (A.+B).^2 = A.^2 .+ 2A.*B .+ B.^2:
t=0:1/200:60;
gauss1 = 250*normpdf(t,19,0.84);
gauss2 = 250*normpdf(t,22,0.87);
m1 = gauss1.^2+2*gauss1.*gauss2+gauss2.^2;
m2 = (gauss1+gauss2).^2;
fprintf('Difference: %1.2f\n', sum(abs(m1-m2)));
subplot(2,1,1);plot(m1);title('Method 1');
subplot(2,1,2);plot(m2);title('Method 2');

I get:
Difference: 0.00

and

So, the math is correct. You may still be wondering, how come the third Gaussian $G_1G_2$ is not apparent in the plot? This is the point where it's useful to challenge your assumption that the third Gaussian should be visible. If we plot $G_1G_2$ by itself we obtain:
plot(2*gauss1.*gauss2);


Look at the amplitude and the position in the x-axis: this third Gaussian is not apparent in the plot because it is completely swamped by the other two, much larger, surrounding peaks.