Error using bvp4c: "Warning: Unable to meet the tolerance without using more than 5000 mesh points. The last mesh of 4921 points and the solution are available in the output argument."

Hello,

I'm using bvp4c to solve for a reaction-diffusion equation. When I set the initial concentration to 100, I get this error saying

"Warning: Unable to meet the tolerance without using more than 5000 mesh points.
 The last mesh of 4921 points and the solution are available in the output argument.
 The maximum residual is 0.512336, while requested accuracy is 0.001. 
> In bvp4c (line 323)
  In steady_state_bvp4c_FINAL (line 38) 
Warning: Imaginary parts of complex X and/or Y arguments ignored 
> In steady_state_bvp4c_FINAL (line 44) 
Warning: Imaginary parts of complex X and/or Y arguments ignored 
> In steady_state_bvp4c_FINAL (line 50) "

Am I setting my boundary conditions incorrectly? I want the left hand side (at x = 0, C_L = C_L0) and on the right hand side I want the flux to equal to 0 (at x = x_f, dC_L/dx = 0).

function steady_state_bvp4c_FINAL
clear all; clc;

%Diffusion coefficient
D_ij = 30; %3.0*10^-7 cm^2/s -> 30 um^2/s

%Initial concentration at x=0
L0 = 100;

%Total length
x_f = 3500;

%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;

L = L0./2;
solinit = bvpinit(linspace(0,x_f,11),[L 0]);

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);

figure(1)
plot(sol.x,(sol.y(1,:)),'LineWidth',1)
xlabel('Distance (\mum)')
ylabel('Concentration (nM)')
axis([0 x_f 0 L0])

figure(2)
plot(sol.x,(sol.y(1,:))./L0,'LineWidth',1)
xlabel('Distance (\mum)')
ylabel('Normalized Concentration (C_L/C_L_0)')
axis([0 x_f 0 1])

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_1) + sqrt((-1 - (C_L./K_1)).^2 -4.*(-(2./K_4) - ((2.*C_L)./(K_2.*K_4)) - ((2.*C_L)./(K_2.*K_4.*K_i)) - ((2.*C_L.^2)/(K_2.*K_3.*K_4.*K_i))).*N_T))/(4.*((1./K_4) + (C_L./(K_2.*K_4)) + (C_L./(K_2.*K_4.*K_i)) + (C_L.^2./(K_2.*K_3.*K_4.*K_i))));

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 = N_T-(2.*N_R2)-(4.*N_pR)-N_R-(4.*N_rR)-(2.*N_r2);
R = N_T-(2.*N_R2)-(4.*N_pR)-N_r-(4.*N_rR)-(2.*N_r2);
r_2 = (N_T./2)-(N_R./2)-N_R2-(2.*N_pR)-(N_r./2)-(2.*N_rR);
rR = (N_T./4)-(N_R./4)-(N_R2./2)-N_pR-(N_r./4)-(N_r2./2);
pR = (N_T./4)-(N_R./4)-(N_R2./2)-(N_r./4)-N_rR-(N_r2./2);
R_2 = (N_T./2)-(N_R./2)-(2.*N_pR)-(N_r./2)-(2.*N_rR)-N_r2;

R_L = (2.*d_1.*(R))-(2.*k_1.*C_L.*r)+...
    ((2.*d_2.*(rR)))-((k_2.*C_L.*(r_2)))+...
    (d_3.*(R_2))-(2.*k_3.*C_L.*(pR));

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

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

Answer 1

[Torsten]

If the code works for values of L0 smaller than 100, try approaching 100 by subsequently calling bvp4c with values L0start, L01, L02,...,L0=100.

Take a look at

https://de.mathworks.com/matlabcentral/answers/428271-how-to-write-matlab-bvp4c-code-for-moving-boundary-surface-problem

to see how this can be implemented.

Best wishes

Torsten.