% Calculate the potential energy of a unit charge distributed as a 
% series of point charges along the x axis, with positions given
% by array xs. It is assumed that xs has at least 2 elements. 
function U = energy(xs)

N = numel(xs); % Number of charges
Qtot = 1.0; % 1 unit of charge for the whole rod.
Q = Qtot / N; % charge of an individual point charge.

U=0;           % Initialize energy at zero
for i=1:N-1
   U = U + Q^2 * sum(1./abs(xs(i) - xs(i+1:end)));
end
