After reviewing your code, I was able to figure out what was troubling you. The flow of your code and equations are all correct, however, you're making a small mistake when creating the normal distribution function (fun = @(x) G;).
Basically what is happening is that you have defined the variable x as a vector of random number already. These vector values are incorrectly used in p1, where x is supposed to be a function variable. As a result, G is just a constant vector. The output of a definite integral should be a scalar value in this case (around 68% as you mentioned) and not a vector.
Fortunately, the fix is quite straightforward, you just need to remove the definition of x as f and then fix your function. Here is how I would implement the integral:
% Generate normally distributed random numbers
nSamples = 100000;
f = 1*randn(nSamples,1);
% Compute normal distribution statistics
m = mean(f);
s = std(f);
% Define normal distribution and integrate over [-1,1] interval
normFun = @(x) 1/(sqrt(2*pi)*s)*exp(-(x-m).^2./(2*s^2));
cdf = integral(normFun,-1,1)
Hope this helps to clarify what was happening!
