Help with a complex 2nd order, nonlinear ODE BVP with complex boundary conditions

I am attempting to solve the Poisson-Boltzmannn equation, which describes the distribution of potential and concentration in relation to a charged surface in a liquid, for a cylinder. The equation this presents is a 2nd order nonlinear ODE that cannot be represented as an IVP, but instead must be solved as a BVP. The initial ODE from this relationship is:

y''=2*z*F*C_bulk/ereo/Vt*exp(mu_bulk-mu_x)*sinh(z*yd)-y'/x

I know in order to use bvp4c in Matlab I have to break this into a system of first order ODEs, so I have set y1=y' to generate the following system:

(1) y1=y'

(2) y1'=y2=2*z*F*c_bulk/ereo/Vt*exp(mu_bulk-mu_x)*sinh(z*y_d)-y1/x

with the boundary cnoditions:

(1) y1 @ x(1,end) = 0

(2) y @ x(1,1) = y_d

There are additional equations that define the system, which include:

mu_x = eta_x.*(8-9*eta_x+3*eta_x.^2)./((1-eta_x).^3)

eta_x = vion*(c_cation+c_anion)

c_cation = c_bulk*exp(-z*y+mu_bulk-mu_x)

c_anion = c_bulk*exp(z*y+mu_bulk-mu_x)

c_x = c_cation+c_anion

q_x = z*c_cation-z*c_anion

c_tot = sum(c_x*pi*(x(1,n+1).^2-x(1,n).^2))

q_tot = sum(q_x*pi*(x(1,n+1).^2-x(1,n).^2))

y_d = Vcell/2/Vt-q_tot*dion*F/2/ereo/Vt

In these equations, variables directly related to distance from the charged surface, x, are row vectors with identical length to the vector, x. These variables include: y1, y2, c_cation, c_anion, c_x, q_x, mu_x, and eta_x.

Some values are constant and they are included in the code I have already generated which will be posted below.

I have generated a code that will predict reasonable values for y if given reasonable mu_x and y_d values (see the original system of ODEs), but the issues I am having are:

(1) How do I include the additional equations in the existing model,and

(2) These equations must be solved iteratively and both depend on and affect the solution of the system of ODEs. How do I perform this iterative calculation?

Here is the code I have generated:

-----------------------------------------------------------------------------------------

function dydx = mybvp(x,y)
global F z c_bulk ereo Vt mu_bulk mu_x
dydx = [ y(2,:).*10^-12; 2*z*F*c_bulk/ereo/Vt*exp(mu_bulk-mu_x(1,:))*sinh(y_d)];

    end

function in = mybvpin (x)
    in=[(9*10^-6.*((x*10^12)^2))-(0.0137.*(x*10^12))+(y_d)
       (1.8*10^-5.*(x*10^12))-0.0137];
end

clear; clc;

global MW %make the variable, MW, usable in all functions
MW=60;%input('Enter the molar mass of the salt ion in g/mol '); %retrieve the molar mass from user

global z %make the variable, z, usable in all functions
z=1;%input('Enter the valence number of the salt ions '); %retrieves the valence number of the salt

global c_bulk %make the variable, c_bulk, usable in all functions
c_bulk=500/MW; %input('Enter the bulk salt concentration in mg/L ')/MW; %retrieve the bulk concentration from user and convert to mol/m3

global dion %make the variable, dion, usable in all functions
dion=700*10^-12; %input('Enter the hydrated ionic diameter in pm ')/10^12; %retrieve the hydrated ionic diameter and convert it to m

global rion %make the variable, rion, usable in all functions
rion=dion/2; %calculate the hydrated ionic radius

global del %make the variable, del, usable in all functions
del=rion; %set the Stern layer thickness

global d %make the variable, d, usable in all functions
d=1*10^-12; %input('Enter the pore diameter in nm ')/10^9; %retrieve the pore diameter and store it as d in units of m

global Vcell %make the variable, Vcell, usable in all functions
Vcell=1.2; %input('Enter the applied voltage in volts '); %retrieve the applied voltage

global vion %make the variable, volume, usable in all functions
vion=pi/6*dion^3; %calculate the volume of one ion from its diameter

global Av %make the variable, Av, usable in all functions
Av=6.022*10^23; %set Avogardo's number

global eta_bulk %make the variable, eta_b, usable in all functions
eta_bulk=vion*c_bulk*2*Av; %calculate the volume fraction of ions

global mu_bulk %make the variable, mu_bulk, usable in all functions
mu_bulk=eta_bulk*(8-9*eta_bulk+3*eta_bulk^2)/((1-eta_bulk)^3); %calculate the bulk electrochemical potential correction in accordance with the CS correction

global e %make the variable, e, usable in all functions
e=-1.60218*10^-19; %set the charge of an electron

global ereo %make the variable, ereo, usable in all functions
ereo=78.64*8.854*10^-12; %set the dielectric permittivity of water

global Vt %make the variable, Vt, usable in all functions
Vt=0.025692577; %set the thermal voltage

global F %make the variable, F, usable in all functions
F=96485; %set the Faraday costant

global y_0 %make the variable, y_0 usable in all functions
y_0=Vcell/2/Vt; %Input the equation for the potential at the charged surface

solinit = bvpinit(linspace(0,155*10^-12,100),@mybvpin);
sol = bvp4c(@mybvp,@mybvpbc,solinit);
x=linspace(0,155*10^-12,20);
y=deval(sol,x);
plot(x,y(1,:));

Thanks in advance for your assistance!

[Torsten]

Could you clarify where the additional equations depend on the solution variables y1 and y2 ?

Best wishes

Torsten.

[Whitewater]

I'm sorry. There was a typo in my original equations. y_d in the equations for c_cation and c_anion should instead be y, which is also a vector with the same length as x.

c_cation = c_bulk*exp(-z*y+mu_bulk-mu_x)

c_anion = c_bulk*exp(z*y+mu_bulk-mu_x)

Also, perhaps I am a little confused about what is being returned by the function. Above the equations are dependent on the local potential, y, which I believe is being returned by the bvp4c routine in the first row with how I have the problem set up. This is y where y1=y' and y2=y''.