bvp4c setup for second order ODE
4 views (last 30 days)
Show older comments
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/191332/image.png)
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!
0 Comments
Answers (1)
Torsten
on 18 Jun 2018
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);
9 Comments
Torsten
on 19 Jun 2018
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
on 22 Oct 2020
Wow this was incredibly helpful as I faced a similar problem. Thank you!
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!