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

Hello,

I am working on solving a system of 3 ODEs using the bvp4c solver. I am receiving the common error "Unable to solve the collocation equations -- a singular Jacobian encountered." From what I have read, this is an error that can be due to a number of problems with the the equation function, boundary value function, or the initial guess. I have tried altering a few things but to no avail.

The .pdf attached details the system I am working with as well as the boundary conditions. I am working to recreate a figure from a research paper and have attached that figure which should give an idea what the solution hopefully will look like.

I believe that the error is likely in my initial guess as I am not quite sure on how to best create that initial guess. I have tried to use functions for the guess that roughly match the solutions plotted in the figures I mentioned above.

Any help or insight is greatly appreciated!

**update July 1st: found a sign error in the first boundary condition and flipped the subraction to addition

global d alpha phi D_0 D_I D_I3 k_edr k_eer k_rg k_ext I2_0 I_0 Keq b gamma theta beta m1 n_0 I3_0 e0 

e0 = 1.602e-19; % elementary charge (Coulombs)
d = 12.5e-4;    % cm
alpha = 0.2e-4; % cm^-1
phi = 6.1e16;   % cm^-2 s^-1
D_0 = 0.4;      % cm^2 s^-1
D_I = 4.4e-6;   % cm^2 s^-1
D_I3 = 3.6e-6;  % cm^2 s^-1
k_edr = 3e-11;  % cm^3 s^-1
k_eer = 2e-16;  % cm^3 s^-1
k_rg = 1.3e-15; % cm^3 s^-1
k_ext = 10^12;  % cm s^-1
n_0 = 1e14;     % cm^-3
I2_0 = 0.05;    % M
I_0 = 0.6;      % M
Keq = 10^8.3;
b = 1;           
gamma = 1.4;    % gamma/omega^2  
theta = 0.45;
beta = 0.35;
m1 = 1.55;

% solving for the equilibrium concentrations of I3_0 and I_0
syms I3
eqn = solve(I3/((I_0 - I3)*(I2_0 - I3)) == Keq, I3);
I3_0 = double(eqn(1));
I_0 = I_0 - I3_0;

x = linspace(0,d,51);
solinit = bvpinit(x,@guess);
sol = bvp4c(@bvfcn,@bc,solinit);

%---Functions----------------------------------------
% u(1) = n, u(2) = dn/dx, u(3) = [I-], u(4) = d[I-]/dx, u(5) = [I3-],
% u(6)= d[I3-]/dx
%
% dudx(1) = dn/dx, dudx(2) = d2n/dx2, dudx(3) = d[I-]/dx, 
% dudx(4) = d2[I-]/dx2, dudx(5) = d[I3-]/dx, dudx(6) = d2[I3-]/dx2

function dudx = bvfcn(x,u)
global phi alpha k_eer k_edr k_rg gamma theta D_0 D_I D_I3

S = phi.*alpha.*exp(-alpha.*x)./(k_rg.*u(3)+k_edr.*u(1));
dudx = [u(2)
        (-phi.*alpha.*exp(-alpha.*x)+k_eer.*u(5).*u(1)+k_edr.*S.*u(1))./(gamma.*(1-theta).*D_0)
        u(4)
        (-1.5.*k_eer.*u(5).*u(1) + 1.5.*k_rg.*S.*u(3))./(gamma.*theta.*D_I./(1.5))
        u(6)
        (0.5.*k_eer.*u(5).*u(1) + 0.5.*k_rg.*S.*u(3))./(gamma.*theta.*D_I3./(1.5))];
end

function res = bc(ul,ur)
global e0 theta D_0 k_ext gamma
j = e0.*k_ext.*ul(1);

res = [ul(2)+(j./(e0.*(1-theta).*D_0.*gamma))
      ur(2)
      ul(4)
      ur(4)+(1.5.*j./(e0.*(1-theta).*D_0.*gamma))
      ul(6)
      ur(6)-(0.5.*j./(e0.*(1-theta).*D_0.*gamma))];
end

function g = guess(x)
global n_0 I_0 I3_0

g = [n_0.*sqrt(x.*10^3)
    n_0./sqrt(x.*10^3)
    I_0.*(2+cos(x.*10^3+pi))
    sin(x.*10^3)
    I3_0.*cos(x.*10^3)
    -sin(x.*10^3)];
end 

[Alan Stevens]

There seems to be an inconsistency between the right side boundary conditions and the graph. The slope of I- at x = d is <= 0 (according to the constants specified in your program), yet it is clearly positive on the graph. The opposite is true of I3-. Am I missing something here?

[Wes Fermanich]

Ah yes sorry, I should have clarified. In plot a, d corresponds to the dashed line where x = 12.5 um. So for this instance I am looking just at the part of the graph that is left of the dashed line. It doesn't look readily true that the slope of I- and I3- are equal to zero at those points but I believe that there is supposed to be type of inflection-like point near the dashed line for both lines.