I would do something like this —
x1 = 0:0.5:10;
x2 = x1 + 2;
x3 = x1 - 0.5;
y1 = 1 + 0.2*x1 + randn(size(x1))*0.2;
y2 = 1.5 + 0.1*x2 + randn(size(x2))*0.1;
y3 = 0.5 + 0.3*x3 + randn(size(x3))*0.3;
DM = [x1, x2, x3;
ones(size([x1, x2, x3]))].';
B = DM \ [y1, y2, y3].'
LBF = DM * B;
figure
plot(x1, y1, '.')
hold on
plot(x2, y2, '.')
plot(x3, y3, '.')
plot(DM(:,1), LBF, '-r')
hold off
legend('y_1','y_2','y_3','Line of Best Fit', 'Location','best')
text(1, 3, sprintf('y = %.3f\\cdotx%+.3f',B))
The strategy here is to concatenate the independent variable vectors and the dependent variable vectors separately, then use them as ordinary independent and dependent variable vectors in the regression (linear here for simplicity). The independent variable vectors can be different values and different lengths (although the same lengths here) making this approach reasonably robust.
A nonlinear approach would proceed similarly, however if different parameters were to be estimated for different independnet and dependent variable vectors, all the sizes would have to be the same for all the vectors, requiring a different approach.
.
