The objective funciton could be something like this —
gausfit = @(b,x) b(1) .* exp(-(x-b(2)).^2/b(3));
y = gausfit([2,4,3],x) + randn(size(x))*0.5;
[B,fv] = fminsearch(@(b)norm(y-gausfit(b,x)), rand(3,1))
plot(x, gausfit(B,x), '-r')
Use the peak locations, amplitudes, and widths, and use a for loop to fit each peak. Much depends on how well the peaks are defined in the data, so it may not be possible to fit all of them accurately.