bvp4c setup for second order ODE

Hello, I'm trying to solve for this second order ODE in steady state using bvp4c with the boundary conditions where at x=0, C_L=1 and x=100, C_L=0.

function bvp4c_test
%Diffusion coefficient
D_ij= 1*10^-6;

%Rate constants
k_1 = 0.00193;
k_2 = 0.00255;
k_3 = 4.09;
d_1 = 0.007;
d_2 = 3.95*10^-5;
d_3 = 2.26;

%Equilibrium constants
K_1 = 3.63;
K_2 = 0.0155;
K_3 = 0.553;
K_4 = 9.01;
K_i = 0.139;

%Total Receptors
N_T =1.7;

solinit = bvpinit(linspace(0,10,11),[1 0]);

sol = bvp4c(@odefcn,@twobc,solinit);

plot(sol.x,sol.y(1,:))

function dC_Ldx = odefcn(C_L,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T)
dC_Ldx = zeros(2,1);

dC_Ldx(1) = C_L(2);

N_r = (-(1+(C_L./K_2))+...
    sqrt(((1+(C_L./K_2)).^2)-4.*((2./K_4).*(1+(C_L./K_2).*(1+(1./K_i).*(1+(C_L./K_3))))).*(-N_T)))./...
    (2.*(2./K_4).*(1+(C_L./K_2).*(1+(1./K_i).*(1+(C_L./K_3)))));

N_R = (C_L./K_1).*N_r;
N_r2 = (1./K_4).*(N_r).^2;
N_rR = (C_L./(2.*K_2.*K_4)).*(N_r).^2;
N_pR = (C_L./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
N_R2 = ((C_L.^2)./(K_2.*K_3.*K_4.*K_i)).*(N_r).^2;

R_L = ((2.*d_1.*(N_T-N_r-N_r2-N_rR-N_pR-N_R2))-(2.*k_1.*C_L.*N_r))+...
    ((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*C_L.*(N_T-N_r-N_R-N_rR-N_pR-N_R2)))+...
    (d_3.*(N_T-N_r-N_R-N_r2-N_rR-N_pR))-(2.*k_3.*C_L.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));

dC_Ldx(2) = -R_L./D_ij;
end

function res = twobc(ya,yb)
res = [ ya(1)-1; yb(1) ];
end
end

When I run it, I encounter the error stating Attempted to access C_L(2); index out of bounds because numel(C_L)=1.

| | _Error in bvp4c_test/odefcn dC_Ldx(1) = C_L(2);

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

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

Error in bvp4c_test sol = bvp4c(@odefcn,@twobc,solinit);_||

Where am I going wrong? I am still very new to MATLAB and appreciate any help that I can get. Thank you!

Answer 1

[Torsten]

sol = bvp4c(@(x)odefcn(x,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T),@twobc,solinit);

[Jennifer Yang]

With that, I get the errors

Error using bvp4c_test>@(x)odefcn(x,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T) Too many input arguments.

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

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

Error in bvp4c_test sol = bvp4c(@(x)odefcn(x,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T),@twobc,solinit);

[Torsten]

Sorry, should be

sol = bvp4c(@(x,y)odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T),@twobc,solinit);

and

function dC_Ldx = odefcn(x,C_L,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T)

[Jennifer Yang]

Thank you for your help so far! Question, why does the sol need to have a y when there isn't y anywhere used in this?

It seems like with each new fix, there is a new error...

Error using bvp4c Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in bvp4c_test sol = bvp4c(@(x,y)odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T),@twobc,solinit);

[Torsten]

R_L is a 2-element vector. Thus the assignment

dC_Ldx(2) = -R_L./D_ij;

does not make sense.

[Jennifer Yang]

So is the problem with the way R_L is written or the way I'm defining my second order differential equation?

[Torsten]

Following your formula for N_r and R_L, you will have to replace C_L in these expressions by C_L(1).

[Jennifer Yang]

Sorry if I am misunderstanding, are you saying replace all the C_L in my N_r and R_L equations with C_L(1) or after the N_r and R_L code, put C_L=C_L(1)?

[Torsten]

function dy = odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T)

C_L    = y(1);
dC_Ldx = y(2);

N_r = (-(1+(C_L./K_2))+...
    sqrt(((1+(C_L./K_2)).^2)-4.*((2./K_4).*(1+(C_L./K_2).*(1+(1./K_i).*(1+(C_L./K_3))))).*(-N_T)))./...
    (2.*(2./K_4).*(1+(C_L./K_2).*(1+(1./K_i).*(1+(C_L./K_3)))));

N_R = (C_L./K_1).*N_r;
N_r2 = (1./K_4).*(N_r).^2;
N_rR = (C_L./(2.*K_2.*K_4)).*(N_r).^2;
N_pR = (C_L./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
N_R2 = ((C_L.^2)./(K_2.*K_3.*K_4.*K_i)).*(N_r).^2;

R_L = ((2.*d_1.*(N_T-N_r-N_r2-N_rR-N_pR-N_R2))-(2.*k_1.*C_L.*N_r))+...
    ((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*C_L.*(N_T-N_r-N_R-N_rR-N_pR-N_R2)))+...
    (d_3.*(N_T-N_r-N_R-N_r2-N_rR-N_pR))-(2.*k_3.*C_L.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));

dy = zeros(2,1);
dy(1) = dC_Ldx;
dy(2) = -R_L/D_ij;

end

[Utkarsh Tiwari]

Wow this was incredibly helpful as I faced a similar problem. Thank you!