Use the r values as an indicator of which distribution f or g to draw from. For those r < 0.4, generate normal random numbers from f, and generate the rest from g.
And since both f and g are normal with the same sigma of 1 but differing only in mu of 1, you could just do this (the rand( ) part simply adds a mu of 1 to those numbers to be drawn from f):
result = randn(1,500) + (rand(1,500) < 0.4);
In general, if you had two normal density functions f and g, and they had sigmas sf and sg and mus mf and mg, and you wanted to combine them with factors pf and pg (pf + pg = 1) and generate n random numbers, it would look something like this:
r = rand(1,n) < pf;
result = randn(1,n);
result(r) = result(r)*sf + mf;
result(~r) = result(~r)*sg + mg;
Or maybe someone out there knows of a function in the Statistics Toolbox that can do this.