Interesting. The 2020b result seems to be correct, based on the formula using sums of squares. (See, e.g. the definition on wikipedia page.)
% Original data
data =[1.8286, 2.0125; 2.5571, 2.3886; 2.0747, 2.3476; 1.6216, 2.0296;
4.6010, 2.4861; 2.6089, 2.4834; 2.0093, 2.5387; 2.4757, 2.5660;
2.5475, 2.6642; 2.4854, 2.6650];
% Define x,y for convenience
x = data(:,1);
y = data(:,2);
mdl = fitlm(x,y,'Intercept',false);
R2 = mdl.Rsquared.Ordinary
% Predicted y value, at x data
y_pred = predict(mdl,data(:,1));
% Calculate R^2 via sums of squares of residuals and total
R_squared_via_sum_of_squares = 1 - sum((y - y_pred).^2)/sum((y - mean(y)).^2)
Based on some playing around, but not quite well formulated enough to post here, I think 2020a might have been using a formula based on the correlation of y data vs. y predicted that is only valid when an intercept term is estimated, and therefore not accurate in your case.
