Just choose any arbitrary vector for ‘z’ that is long enough to completely support the signals to be convolved. For a Gaussian, that is relatively easy to describe. (For other signals, it may be less obvious.)
z = (min(c1,c2)-10):(max(c1,c2)+10);
x = exp(-(z-c1).^2/(2*s1.^2));
h = exp(-(z-c2).^2/(2*s2.^2));
y = my_convolution(x, h);
same_range = (fix(numel(y)/2)-fix(numel(z)/2))+(1:numel(z));
comp = conv(x, h, 'same');
function y = my_convolution(x, h)
y(i) = y(i) + X(j)*H(i - j + 1);
The ‘same_range’ index vector appears to work reasonably well.
Experiment to get different results.