function [sol] = pbeexp()
el = 4.8*10^-10;
eps = 80;
k = 1.3807*10^-16;
T = 293 ;
l = el^2/(eps*k*T);
pn = k*T/el;
psiwall =300*10^-3*0.00333564;
x=linspace(0,35*10^-7,10000)/l;
solinit = bvpinit(x,[psiwall/pn 0]);
options = bvpset('RelTol',1e-5,'AbsTol',1e-5);
sol = bvp4c(@pbeode,@pbebc,solinit,options);
xint=linspace(sol.x(1),sol.x(end),size(sol.x,2));
sxint=deval(sol,xint);
dist = xint;
pot = sxint(1,:);
plot(dist,pot)
function dydx = pbeode(x,psi)
el = 4.8*10^-10;
eps = 80;
k = 1.3807*10^-16;
T = 293 ;
l = el^2/(eps*k*T);
H = 7.8*10^-7;
a3 = (1/4.6)*10^3 ;
czero = 6.023*10^20*0.1 ;
nu = 2*czero*a3;
pn = k*T/el;
A1 = (4*pi*czero*l*sinh(psi(1)))/(1+2*nu*((sinh(psi(1)/2)^2)));
dydx = [ psi(2)
(2*l^2*A1)];
end
function bc = pbebc(psia,psib)
psiwall = 300*10^-3*0.00333564;
pn = k*T/el;
bc = [ psia(1)-psiwall/pn
psib(1)];
end
end