function y=gafit(x0) % this function takes in one initial guess
                     % vector of length 2 x0 = <x0(1), x0(2)>
                     % corresponding to the unknowns x = <A, B>

load gaussfit.dat % data to be fit by a gaussian
x=gaussfit(:,1);  
y=gaussfit(:,2);

% Minimizing the following sum is equivalent to minimizing
% the mean root square error
sum=0;
for j=1:length(x)
   sum=sum+ ( x0(1)*exp(-x0(2)*x(j)^2) ...
      -y(j)  )^2; 
end
y=sum;

% or instead of the 'for loop' this could be faster:
%y=sum( ( x0(1)*exp(-x0(2)*x.^2)-y ).^2 );