The Statistics Toolbox has a function lsline that does this automagically but for some reason that is totally unexplicable to me, they've neutered it so it ignores any/all lines that are not just points/markers. Dunno why they think use cases like this aren't common as well.
But a workaround could be
hAx=gca; hL=hAx.Children; % retrieve the axes and line handles
hF=arrayfun(@(h)fitlm(h.XData,h.YData,"linear"),hL,'Uni',0); % fit each in turn, save fit handles
LS=get(hL,'Linestyle'); % get the used linestyle to put back
MK=get(hL,'Marker'); % and marker, too
set(hL,{'Linestyle'},{'none'},'Marker',{'o'}) % turn line off, add marker
hLSqL=lsline(hAx); % now add the LS line
set(hL,{'Linestyle'},LS,{'Marker'},MK) % put back the original line, marker
% optional instead to keep a marker -- can be hard to see with both lines where data are
%set(hL,{'Linestyle'},LS,{'Marker'},{'none'}) % put back the original line, leave marker
Alternatively, you can (and I'd recommend) use the X,Y data you had when you drew the plot and follow a similar route but use those data directly.
