1 missing bracket and you missed to define dx:
d = 30;
E0=8.85e-12;
c=.1;
dx = 0.1;
a = 0; b = 1;
x = a:dx:b;
N = length(x);
f = @(x) ((d)./(4.*pi.*E0))*(c./((sqrt(c.^(2)+x.^(2))).^(3)));
y = f(x);
S = 0;
for k = 1:N-1
S = S + 0.5*(y(k+1) + y(k))*(x(k+1) - x(k));
end
S_trap = trapz(x,y);
S_int = integral(f,a,b);
% print results to pdf
fprintf('INTEGRAL #1\n')
fprintf('Matlab''s integral function yields: %.3f\n',S_int);
fprintf('Matlab''s trapezoid function yields: %.3f\n',S_trap);
fprintf('My trapezoid method yields: %.3f\n\n',S);
% plot integrand and print my trapezoid integration value
figure(1); fplot(f,[a,b])
xlabel('x'); ylabel('f(x)')
title([{'Coulombs Law'},...
'The approximate area using my trapezoid method is: ',...
num2str(S,6)]);
grid on
