Issue with PDE Non-constant Boundary Condition

I am using solvepde to obtain a solution to the PDE where It has a Dirichlet boundary condition on the left, a section of non-constant Neumann boundary condition on the right and zero flux everywhere else. The non-constant boundary condition is defined by where was previously solved for and all other variables are constants. I have defined this non-constant boundary condition with bcfun in the code below. Using the default initial condition, I get the error message "Error using pde.EquationModel/solveStationaryNonlinear (line 27) Unsuitable initial guess U0 (default: U0=0)." Since is nearly constant with y, (it is only different at the ends where the boundary condition changes from non-constant to zero flux) I got approximate results by assuming it was constant, and I set an initial condition based off these results. However, the two figures show that the results of the code below are very different from the results of the approximation, which I would not expect since is almost constant everywhere on the boundary.

Does anyone know what the problem is or how to fix it?

Here is the full code for reference:

% Define constants
e = 1.60217662*10^-19;
l = 10*10^-6;
sige = 3.37*10^-4;
sigi = 18;
T = 1000;
R = 8.314;
kb = 1.381*10^-23;
h = 6.626*10^-34;
F = 96485;
n = -0.02;
c = 1;
pH2 = 0.2;
pH2O = 1-pH2;

% Determine value of muO2 at the right boundary
G = -247500+55.85*T;
keq = exp(-G/(R*T));
pO2 = (pH2O/pH2)^2*(1/keq^2);
mu2 = kb*T*log(pO2);

% Rectangle is code 3, 4 sides, followed by x-coordinates and then y-coordinates
R1 = [3,4,0,l,l,0,0,0,.025,.025]';
R2 = [3,4,0,l,l,0,0.025,0.025,.075,.075]';
R3 = [3,4,0,l,l,0,0.075,0.075,.1,.1]';
geom = [R1,R2,R3];
% Names for the two geometric objects
ns = (char('R1','R2','R3'))';
% Set formula
sf = 'R1+R2+R3';
% Create geometry
g = decsg(geom,sf,ns);

% Create mu geometry model
mumodel = createpde;
geometryFromEdges(mumodel,g);

% Apply boundary conditions
applyBoundaryCondition(mumodel,'dirichlet','Edge',8,'u',0);
applyBoundaryCondition(mumodel,'dirichlet','Edge',9,'u',0);
applyBoundaryCondition(mumodel,'dirichlet','Edge',10,'u',0);
applyBoundaryCondition(mumodel,'dirichlet','Edge',6,'u',mu2);
applyBoundaryCondition(mumodel,'neumann','Edge',1,'g',0);
applyBoundaryCondition(mumodel,'neumann','Edge',3,'g',0);
applyBoundaryCondition(mumodel,'neumann','Edge',4,'g',0);
applyBoundaryCondition(mumodel,'neumann','Edge',2,'g',0);
applyBoundaryCondition(mumodel,'neumann','Edge',5,'g',0);
applyBoundaryCondition(mumodel,'neumann','Edge',7,'g',0);

% Solve PDE for mu and plot results
specifyCoefficients(mumodel,'m',0,'d',0,'c',1,'a',0,'f',0);
generateMesh(mumodel,'Hmax',l);
solmu = solvepde(mumodel);
mu = solmu.NodalSolution;
f1=figure;
pdeplot(mumodel,'XYData',mu)
title('Mu vs. x and y')
xlabel('x-position')
ylabel('y-position')

% Create phi geometry model
phimodel = createpde;
geometryFromEdges(phimodel,g);

% For the approximation I defined the constant Neuman boundary condition:
% A = (sigi/(4*e)*evaluateGradient(solmu,l,0.05)+2*F*c*pO2^n)/sigi;
% and substituted A for bcfun in the boundary condition applications below

% Define nonconstant Neumann boundary condition
Ie = -2*F*(c*pO2^n);
bcfun = @(location,state)(sigi/(4*e)*evaluateGradient(solmu,l,location.y)-Ie)/sigi;

% Make sure initial condition is suitable
setInitialConditions(phimodel,0);
setInitialConditions(phimodel,-0.7,'Edge',6);

% Apply boundary conditions
applyBoundaryCondition(phimodel,'dirichlet','Edge',8,'u',0);
applyBoundaryCondition(phimodel,'dirichlet','Edge',9,'u',0);
applyBoundaryCondition(phimodel,'dirichlet','Edge',10,'u',0);
applyBoundaryCondition(phimodel,'neumann','Edge',6,'g',bcfun);
applyBoundaryCondition(phimodel,'neumann','Edge',1,'g',0);
applyBoundaryCondition(phimodel,'neumann','Edge',3,'g',0);
applyBoundaryCondition(phimodel,'neumann','Edge',4,'g',0);
applyBoundaryCondition(phimodel,'neumann','Edge',2,'g',0);
applyBoundaryCondition(phimodel,'neumann','Edge',5,'g',0);
applyBoundaryCondition(phimodel,'neumann','Edge',7,'g',0);

% Solve PDE for phi and plot results
specifyCoefficients(phimodel,'m',0,'d',0,'c',1,'a',0,'f',0);
generateMesh(phimodel,'Hmax',l);
solphi = solvepde(phimodel);
phi = solphi.NodalSolution;
f2=figure;
pdeplot(phimodel,'XYData',phi)
title('Phi vs. x and y')
xlabel('x-position')
ylabel('y-position')