% Estimate a continuous charge density distribution for a unit 
% charge distributed as a set of N point charges on the x-axis.
function lambda = charge_density(xs)

Npts = numel(xs); % number of point charges
Nint = Npts - 1; % number of intervals between the point charges
dx = 1/Nint; % mean interval size
lambda0 = 1.0; % mean charge density.

% Approach: Divide the mean charge density by a stretching factor
% appropriate to each interval.

% Use a centered difference to estimate stretching for interior points:
lambda(2:Npts-1) = lambda0 * 2*dx ./ ( xs(3:Npts) - xs(1:Npts-2) );
% (truncation at second order)

% Stretching at ends cannot be done with a centered difference:
lambda(1) = lambda0 * dx / ( xs(2) - xs(1) );
lambda(Npts) = lambda0 * dx / ( xs(Npts) - xs(Npts-1) );
% (truncation at first order --- more error prone)
