Simultaneous fitting of multiple datasets?

I'm trying to fit multiple datasets which have the same underlying behavior - there are several gaussian peaks, and the center/width of each peak is the same in all of the datasets. However, the peak height can vary between datasets.

I'm interested in finding the peak centers. Because the fit function, centers, and widths are the same for all datasets, I'd like to fit multiple datasets at the same time, using peak centers as a global variable.

Below is some code which simulates my situation: I end up with three datasets (y_1, y_2, y_3), and while I can fit any one of them, I can't see how to fit them all at the same time.

Any help is appreciated.

****************************************

function multifit

% independent variable
xdat=[(0:.5:10)'];

% === model function ===
function y_out=twogauss(x_array,p_in)
  y_out= gauss(x_array,p_in(1:3)) + gauss(x_array,p_in(4:6)) + p_in(7);
end
% ======================

% signal (if it were perfect) - note peaks not the same height in each dataset
cen_a=3.14; wid_a=1.57;
cen_b=7.20; wid_b=1.42;
y_1=twogauss(xdat, [30,cen_a,wid_a,  20,cen_b,wid_b, 1]);
y_2=twogauss(xdat, [20,cen_a,wid_a,  10,cen_b,wid_b, 2]);
y_3=twogauss(xdat, [15,cen_a,wid_a,  30,cen_b,wid_b, 2]);
y=[y_3];  % y_1, y_2, or y_3

% Add some noise to signal
ydat = y + 0.25*((sqrt(y)).*randn(size(xdat)));
edat=.25*sqrt(y);

% fit the data
p_guess=[20,3,1.6,  20,7,1.4,  0];
p_free=[1,1,0, 1,1,0, 1];  % peak widths are fixed
[y_fit,p_out]=fit(xdat,ydat,edat,@twogauss,p_guess,p_free);

% Display fitted coefficients
disp([' Centers : ',num2str([p_out.p(2) p_out.p(5)])])

% plot data/fit (not required)
hold off
errorbar(xdat,ydat,edat,'o');
hold;
plot(xdat,y_fit,'r-');
axis([-2 12 0 40]);

end

****************************