I am trying to solve a boundary-value problem - a system of 3 ODEs - using the bvp4c command. I am returned the warning below regarding the tolerance. The warning states that it needs to use 1666 mesh points to meet my specified tolerance (1E-9). However, as the warning states itself, the simulation consists of more mesh points (2004 mesh points) than the required 1666. I also tried the simulation with many more mesh points (20005 mesh points), and I still got a similar warning message. Why, then, is it not able to meet the tolerance? I included the exact warning message, and my code (in case anyone wanted to try this for themselves - takes about 20sec to run on my laptop), below.
Warning: Unable to meet the tolerance without using more than 1666 mesh points. The last mesh of 2004 points and the solution are available in the output argument. The maximum residual is 0.305596, while requested accuracy is 1e-09. > In bvp4c at 326 In test_model at 123
% ---------------------- Initialization ---------------------------- % Define mesh in each region xinit_region1 = -x_m:h:0; xinit_region2 = 0:h:1; xinit_region3 = 1:h:x_s;
% Since a region may not be divisble by h, ensure that the boundaries % of the region are defined if( xinit_region1(end) ~= 0 ) xinit_region1 = [ xinit_region1 0 ]; end if( xinit_region2(end) ~= 1 ) xinit_region2 = [ xinit_region2 1 ]; end if( xinit_region3(end) ~= x_s ) xinit_region3 = [ xinit_region3 x_s ]; end
if (x <= 1) % algebraic equations for r, b_er, and b_lr [r,~,exitflag] = fzero(@(r) receptor_equation(r,x,y), [-0.0000001 1000]); % find r numerically (can't find explicit solution) if(exitflag ~= 1) error('Abedsays:rNotFound', 'r was not found. Values: exitflag: %d, y(1): %f, y(3): %f', exitflag, y(1), y(3)); end b_er = (k_ER_p*y(3)*r)/(k_ER_m + a_ER); b_lr = (k_LR_p*y(1)*r)/(k_LR_m + a_LR);
% production of species v_L = v_L_max * heaviside(-x); %ligand production v_E = (v_E_max/(1+(k_E*(b_lr)^n_E)))*(heaviside(x)); %expander production else % sink region r = 0; b_er = 0; b_lr = 0; v_L = 0; v_E = 0; end
% Define boundary conditions function res = bc(YL,YR) res = [ YL(2,1) % dl(-xmin)/dx = 0 YR(1,end) % l(xs) = 0 YL(4,1) % de(-xmin)/dx = 0 YR(3,end) % e(xs) = 0 YL(6,1) % d(b_el(-xmin))/dx = 0 YR(5,end) % b_el = 0 YR(1,1)-YL(1,2) % continuity of l at x=0 YR(2,1)-YL(2,2) % continuity of l' at x=0 YR(3,1)-YL(3,2) % continuity of e at x=0 YR(4,1)-YL(4,2) % continuity at e' at x=0 YR(5,1)-YL(5,2) % continuity of b_el at x=0 YR(6,1)-YL(6,2) % continuity of b_el' at x=0 YR(1,2)-YL(1,3) % continuity of l at x=1 YR(2,2)-YL(2,3) % continuity of l' at x=1 YR(3,2)-YL(3,3) % continuity of e at x=1 YR(4,2)-YL(4,3) % continuity at e' at x=1 YR(5,2)-YL(5,3) % continuity of b_el at x=1 YR(6,2)-YL(6,3) ]; % continuity of b_el' at x=1 end end