This file poorly fits a gaussian (and has a much higher average rating) as compared to others here...
The test main file is wrong: the "rand" noise on the gaussian biases the measurements by 0.5, randn should be used instead.
Same main file used for comparison with another gaussian fit (gaussfit.m):
%% data
x=1:100;
sigma=15; mu=40; A=3;
plot(normpdf(x,mu,sigma))
y=A*exp(-(x-mu).^2/(2*sigma^2))+randn(size(x))*0.5;
hold all;
plot(x,y/(sum(y)),'.');
%% fitting
[sigmaNew,muNew,Anew]=mygaussfit(x,y);
y2=Anew*exp(-(x-muNew).^2/(2*sigmaNew^2));
plot(x,y2/(sum(y2)));
%% bis
[sigmaNew2,muNew2]=gaussfit(x,y/sum(y));
plot(normpdf(x,muNew2,sigmaNew2));
legend('truth','measurements','mygaussfit','gaussianfit')

Great idea,
though I wonder when you do the retransformation from second order polynomial (a0 + a1*x + a2*x^2)to logarithmic gaussian (log(A) - x^2/2*sigma^2 + x*mu/2*sigma^2 - mu^2/2*sigma^2)(line 46 in mygaussfit) shouldn't you calculate the mean from p(2) via mu = A1/-A2 or equivalently mu = A1*2*sqrt(-1/2*A2) instead of mu = A1*2*sigma^2 ?
I might be totally wrong or missed sth, just a quick idea.