BVP solver giving singular jacobian error even with analytical jacobians provided

Hey all, I'm getting a singular Jacobian error while trying to use the solver bvp4c:

??? Error using ==> bvp4c at 252
Unable to solve the collocation equations -- a singular Jacobian encountered.

The only thing is that I'm providing the solver with the analytical jacobians! Any ideas as to what's going on? Here's my code:

One thing I noticed was the the solver always messed up while on x=100.5, with the boundary condition at x=101. Not sure what that indicates.

function potential = poisson(xspace, phia, dphidxa, H, Zn)

    %Random constants stolen from main code. Will clean these up if it proves
    %worthwhile
    ep0 = (8.85e-12)/100;             % permittivity in SI [F/cm]
    ep0 = ep0/10^10;                 %Scaling permittivity down by 10^10
    ep  = 12;                   % relative permittivity
    q   = 1.6e-19;              % charge on an electron [C]
    T  = 463;                 % Temperature, [K]
    kk = 1.381e-23;           % Boltzmann's constant [J/K]
    ni = 5e11;             % number of intrinsic carriers
    ni = ni/10^10;             %Scaling # of intrinsic carriers down by 10^10
    lx = length(xspace);

    %defining the x mesh
    xspace = 1:lx;

      solinit = bvpinit(xspace, @potentialguess);
      options = bvpset('Stats','on', 'FJacobian', @potentialjac, 'BCJacobian', @potentialBCjac, 'RelTol',1e-8);

      sol = bvp4c(@potentialode, @potentialbc, solinit, options);
      xint = xspace;
      Sxint = deval(sol,xint);
      potential = Sxint;

      %dydx sets up the system of first order ODEs, which is how matlab
      %processes these types of problems. y(1) is phia, y(2) is dphidxa. And
      %each entry of this vector is the derivative of each. So the derivative
      %of y(1) is y(2). And the derivative of y(2) is the poisson equation.
      function dydx = potentialode(x,y)
          if floor(x) ~= ceil(x)
              points = [floor(x), ceil(x)];
              Hinterp = interp1(points, H(points), x);
              Zninterp = interp1(points, Zn(points), x);
          else
              Hinterp = H(x);
              Zninterp = Zn(x);
          end
          dydx = [y(2)
                  (q/(ep*ep0))*(ni*exp((q*y(1))/(kk*T)) - ni*exp((-q*y(1))/(kk*T)) + Zninterp - Hinterp)];
      end

      %The residual matrix, which sets the boundary conditions. The computer
      %attempts to make each entry 0 within this matrix during its solution algorithm.
      function res = potentialbc(ya, yb)
         res = [ya(1) - phia(1)
                yb(1) - phia(lx)];
      end

      %The guess matrix provides a guess for the future values of the
      %potential and its derivative. In this case, the guess is simply the
      %old values of the potential and its derivative.
      function guess = potentialguess(x)
          guess = [phia(x)
                   dphidxa(x)];
      end

      %This function is the analytical jacobian - greatly speeds up the
      %working of the solver so the solver can skip calculating it
      %numerically
      function dfdy = potentialjac(x, y)
          dfdy = [0                                            1   ((q^2*ni)/(ep*ep0*kk*T))*(exp((q*y(1))/(kk*T)) + exp((-q*y(1))/(kk*T)))  0];
      end

      function [dBCdya, dBCdba] = potentialBCjac(ya, yb)
          dBCdya = [1  0
                    0  0];
          dBCdba = [0  0
                    1  0];
      end
  end