There were some minor mistakes, see below for the code with inline comments to explain the differences.
A = 1.325;
B = 3.7;
C = 5.74;
% number of e/d values, set to 3 as in the question
nD = 3;
% Relative Roughness range
D = linspace(10.^-6,10.^-2,nD);
% Reynolds number range
E = 108:4000;
% initialize a matrix for the friction factor values, as many rows as the
% ralative roughness range and as many columns as the reynolds number range
f = zeros( nD , numel(E) );
% evaluate the Swamee-Jain formula, in a loop, iterate over the relative
% roughness range
for i = 1:nD
f(i,:) = A ./ log( (1/3.7)*D(i)+(C./E.^0.9) ).^2;
% note that here we use ./ and .^ for element wise operations
end
figure
plot(E, f);
xlabel('Reynolds Number')
ylabel('Friction Factor f')
title('Friction Factor against Reynolds Number')
grid on