Boundary value problem with a changing parameter

Hi!

I'm trying to solve a boundary value problem for a set of odes, with a parameter changing with x.

I have the code set up as below.

But the code would only run if I set AbsTol and Reltol to 1E-1. Also, I'm not getting the expected results.

Could someone please help me out?

Really appreciate!

F=96485;                    
R=8.3144621;          
T=25+273.15;           
Z_OH=-1;
Z_H=1;
LA=1E-4;                       
LC=1E-4;                       
d=1E-5;                        
D_H=5.94/10^10;               
D_OH=3.47/10^10;            

%%%%%%%%%%%%%%%%%%%%
options=bvpset('AbsTol',1e-1,'Reltol',1e-1);
y0=[0.021;258.1;226.0552;486.6069;0;0];
solinit=bvpinit([0,LA+LC+d],y0);
sol =bvp4c(@odefun,@odebc,solinit,options);
%%%%%%%%%%%%%%%%%%%%

p=(linspace(0,LA+LC+d))';
y=(deval(p,sol))';

function [dy]= odefun(x,y)
global F R T  Z_OH Z_H  D_OH  D_H  LA LC d
% y(1)=Electric potential,y(2)=Electric field, y(3)=OH-concentration
% y(4)=Na+ concentration, y(5)=dy(3)/dx, y(6)=dy(4)/dx

if x>=0 && x<=(LA-(5E-6))
    FixedCharge=2000;
elseif x>(LA-(5E-6)) && x<=(LA+(5E-6))
    FixedCharge=1000*(tanh(x*5E5)-tanh((x-LA)*5E5));
elseif x>=(LA+(5E-6)) && x<=(LA+d-(5E-6))
    FixedCharge=0;
elseif x>=(LA+d-(5E-6)) && x<(LA+d+(5E-6))
    FixedCharge=-1000*(tanh((x-LA-d)*5E5)-tanh((x-LA-LC-d)*5E5));
else
    FixedCharge=-2000;
end

Kw=1E-8;

dy = zeros(6,1);

dy(1)=-y(2);
dy(3)=y(5);
dy(4)=y(6);

dy(2)=(D_H*Kw/(D_OH*(y(3)^3))*(2*(y(5)^2)-y(3)*dy(5))+D_H/D_OH*Z_H*F/R/T*Kw/(y(3)^2)*y(5)*y(2)-dy(5)+Z_OH*F/R/T*y(5)*y(2))/(D_H/D_OH*Z_H*F/R/T*Kw/y(3)-Z_OH*F/R/T*y(3));
dy(5)=(dy(6)+Kw*2*(y(5)^2)/(y(3)^3)-Z_OH*F/R/T*((y(6)-y(5)-Kw/(y(3)^2)*y(5))*y(2)+(y(4)+Kw/y(3)-y(3)+FixedCharge)*dy(2)))/(1+Kw/(y(3)^2));
dy(6)=Z_H*F/R/T*(y(6)*y(2)+y(4)*dy(2));
end

function [bc]= odebc(ya,yb)    % Boundary condition
bc=[ya(1)-0.021
    yb(1)-1
    ya(3)-226.0552
    yb(3)-4.4237E-11
    ya(4)-486.6069
    yb(4)-2260.6];
end

[Torsten]

Use the capability of BVP4C to solve multipoint boundary value problems:

https://de.mathworks.com/help/matlab/ref/bvp4c.html#bt5uooc-23

This will cope with the discontinuities of your ODEs because of the FixedCharge parameter.

Best wishes

Torsten.

[Luka]

Hi Torsten, Thank you very much for the help.

I read the link you sent, and it says

"In the boundary conditions function

bcfun(yleft,yright)

yleft(:,k) is the solution at the left boundary of [ak−1,ak]. Similarly, yright(:,k) is the solution at the right boundary of region k. "

But for my case, I don't have boundary conditions for the intervals.

[Luka]

Hi Torsten,

Could you please give me some more advice? I'm still stuck with it.

In the odefun function part, did I express the equations right? Is it ok to have dy/dx at both sides of the equation?

Thank you very much.