function  pb()


%% Applied voltage 
Voltage = 2;       

%% Grid 
Xmax = 10;
Xmin = Xmax/10000;


options = bvpset('RelTol',1e-5,'AbsTol',1e-5,'Nmax',5000);
x = logspace(log10(Xmin),log10(Xmax),100);
x = [0 x];

solinit = bvpinit(x,@pbguess,[],Voltage);


sol = bvp4c(@pbeqns,@pbbc,solinit,options,Voltage);

subplot(2,1,1)
hold off
plot(sol.x,sol.y(1,:),'r');
hold on
plot(sol.x,sol.y(2,:));
xlabel('x')
ylabel('C^+, C_-')

subplot(2,1,2)
plot(sol.x,sol.y(3,:))
xlabel('x')
ylabel('\phi')

%% compute the total charge
charge = trapz(sol.x,(sol.y(1,:)-sol.y(2,:)))/2;

fprintf('V %f q %e q theoretical %e\n',Voltage, ...
	charge,-2*sinh(Voltage/2));




%%%%%%%%%%%%%%%%%%%%%
%% Initial Guess
%%%%%%%%%%%%%%%%%%%%%
function Yinit = pbguess(x,Voltage)

%% The guess is the solution to the linearized, low voltage equations

Phi = Voltage*exp(-x);
Cp  = exp(-Phi);
Cm  = exp( Phi);
E   = -Phi;

Yinit = [Cp Cm Phi E];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Poisson Boltzmann equations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dydx = pbeqns(x,y,Voltage)


%% unpack the incoming vector
Cp   = y(1);
Cm   = y(2);
Phi  = y(3);
E    = y(4);

%% governing equations
dCdx   = [-Cp*E;Cm*E];
dPhidx = E;
dEdx   = (Cm-Cp)/2;


%% repack 
dydx = [dCdx; dPhidx; dEdx];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function res = pbbc(ya,yb,Voltage);

%% unpack the incoming vector
Cpa   = ya(1);
Cma   = ya(2);
Phia  = ya(3);
Ea    = ya(4);

Cpb   = yb(1);
Cmb   = yb(2);
Phib  = yb(3);
Eb    = yb(4);

%% set the BC's
res = [Phia-Voltage;
       Phib;
       Cpb-1;
       Cmb-1];

