N = 2^8 * 100;
dt = 0.00001;

x1 = rand(N, 1);
x2 = rand(N, 1);
q = rand(N, 1) + 1;
v1 = zeros(N, 1);
v2 = v1;

%Start of a timestep
tic
r1 = bsxfun(@minus, x1, x1');
r2 = bsxfun(@minus, x2, x2');
rmag = r1.^2 + r2.^2;
rmag = (rmag.*sqrt(rmag))+1e-4;

r1 = r1./rmag;
r2 = r2./rmag;
r1(1:N+1:end) = 0;
r2(1:N+1:end) = 0;

v1 = v1 - r1*q*dt;
v2 = v2 - r2*q*dt;
x1 = x1 + v1*dt;
x2 = x2 + v2*dt;
toc