% What is the self-energy of a thin conducting rod of unit length and unit charge?
% What is the equilibrium distribution of charge?

% Improved version with equations of constraint.

xmin = -0.5; % coordinate of left end or rod
xmax = +0.5; % coordinate of right end of rod
x0 = (xmin:0.005:xmax)'; % initial charge position array. 
   % Note that sqp uses column vectors when calling @constraints.

% Print out some initial stats...
N = numel(x0)  % Number of charges
Einitial = energy(x0)

maxiter = 1000; % maximum number of iterations
tol = sqrt(eps); % tolerance for convergence

t0 = cputime();

[xfinal, Efinal, info, iter] = sqp(x0, @energy, [], @constraints1, xmin, xmax, maxiter, tol);

execution_time = cputime() - t0

% Print out some final stats...
Efinal
info
iter
fractional_energy_change = (Efinal - Einitial)/Efinal

hold off
plot(x0, xfinal,'r-');
hold on
xlabel('initial x')
ylabel('final x')
title('Test 2 Solution')
print('test2.pdf','-dpdf')
hold off