Help regarding use of BVP4C in solving an ODE.

This is a script I have written :-

clear all
e = 1.602*10^-19; %electron charge
eps = 8.85*10^-12; %permittivity of free space
k = 1.38*10^-23; %boltzmann constant
T = 293 ;%absolute temperature
a3 = 4.6 ; %steric packing
czero = 0.1 ; %bulk counterion concentration
nu = 2*czero*a3;%steric size parameter
l = e^2/(4*pi*eps*k*T);%Bjerrum length
H = 3.8*10^-9;
psival = -10;

newpbe1(l,H,psival)

*Firstly a few constants are defined. Then in the end a function "newpbe1" is called.

The code for the function "newpbe1" is given below. It is saved in a different file. Also the two different functions "pbeode " and "pbebc" are written in two different files and are called in the last line of "newpbe1" where the "bvp4c" function is used.*

function [x,psi] = newpbe1(l,H,psival)

x = [0,logspace(-15,log10(H/2),2e3)]*1/l; % Exponentially distributed grid
solinit = bvpinit(x,[psival 0]); % Initial mesh
options = bvpset('RelTol',1e-7,'AbsTol',1e-7); % Tolerance values

psi = bvp4c(@pbeode,@pbebc,solinit,options);

end

function dydx = pbeode(x,psi,l,czero,nu)
dydx = [ psi(2);(2*l^2*((czero*sinh(psi(1)))/(1+2*nu*(sinh(psi(1)/2))^2)))];
end

function bc = pbebc(psia,psib,psival)
bc = [ psia(1)-psival ; psib(1)];
end

• When I am running this code I am getting the following error message. *

Error using pbeode (line 2)
Not enough input arguments.

Error in bvparguments (line 106)
    testODE = ode(x1,y1,odeExtras{:});

Error in bvp4c (line 130)
[n,npar,nregions,atol,rtol,Nmax,xyVectorized,printstats] = ...

Error in newpbe1 (line 7)
psi = bvp4c(@pbeode,@pbebc,solinit,options);

Error in newpbe (line 14)
newpbe1(l,H,psival)

Kindly help me in resolving the problem. I am unable to fix it. Thanks.

Answer 1

[Mischa Kim]

AGNIVO, check out:

 function my_pde()
 e = 1.602*10^-19; %electron charge
 [...]                        %%remove code for better readability
 psival = -10;
 newpbe1(l,H,psival,czero,nu) %%need to pass all parameter
 end

 function [x,psi] = newpbe1(l,H,psival,czero,nu)
 x       = [0,logspace(-15,log10(H/2),2e3)]*1/l; % Exponentially distributed grid
 solinit = bvpinit(x,[psival 0]); % Initial mesh
 options = bvpset('RelTol',1e-7,'AbsTol',1e-7); % Tolerance values
 par_vec = [l; czero; nu];    %%package parameter in vector
 psi     = bvp4c(@pbeode,@pbebc,solinit,options,par_vec);
 end

 function dydx = pbeode(x,psi,par_vec)
 l     = par_vec(1);
 czero = par_vec(2);
 nu    = par_vec(3);
 dydx  = [ psi(2);(2*l^2*((czero*sinh(psi(1)))/(1+2*nu*(sinh(psi(1)/2))^2)))];
 end

 function bc = pbebc(psia,psib,psival)
 bc = [psia(1) - psival(1); psib(1)]; %%this needs to be a 2-by-1
 end

• I put all code into one MATLAB file (named my_pde.m)
• You need to pass all required parameters to relevant functions
• The output vector in pbebc needs to be a 2-by-1, I simply "fixed" it as shown above. Please verify.

[AG4102M]

Thank you. I have also found the correct solution by myself. But thanks anyway.