Trouble with anonymous function coefficients in PDE Toolbox

I'm running into an issue when I try to use a functional coefficient in my PDE Toolbox problem. In my programmatic script, solvepde returns a solution that is zero everywhere. Curiously, when I enter the same information into the PDE Modeler GUI, I get the solution I expect.

My simulation is a solenoid in axisymmetric cylindrical coordinates; I have a region defining the coil and a region defining the surrounding vacuum. For this problem, the solution u gives the azimuthal component of the magnetic vector potential, m and d are zero, c is 1/mu0, f is the current density and is 1 in the coil and 0 elsewhere, and a is 1/(mu0 r^2). For the MATLAB implementation, I am using x for z and y for r.

Here is my code:

% generate the PDE model
model = createpde;

% specify geometry (in m)
coil_L = 0.04;
coil_IR = 0.01;
coil_thickness = 0.003;
z_max = 0.05;
r_max = 0.03;

% specify material constants and current density
current_density = 1; % A/m

% constants
mu0 = 4*pi*1e-7;

% derived quantities
coil_OR = coil_IR + coil_thickness;

% define the geometry
coil = [3 4 ...
      -coil_L/2   -coil_L/2   coil_L/2    coil_L/2 ...
      coil_IR     coil_OR     coil_OR     coil_IR]';
bounds = [3 4 ...
      -z_max  -z_max  z_max   z_max ...
      0       r_max   r_max   0]';

% combine the geometry elements
geom = [coil, bounds];

% create the geometry
geometryFromEdges(model, decsg(geom));

% apply boundary conditions
applyBoundaryCondition(model,'dirichlet','Edge',[3:5,8],'h',1,'r',0);

% specify coefficients
a_coeff = @(region,state) region.y.^(-2)/mu0;
%a_coeff = @(region,state) zeros(size(region.y));
specifyCoefficients(model, 'm',0,'d',0, 'c', 1/mu0, 'a', a_coeff, 'f', 1, 'Face', 1);
specifyCoefficients(model, 'm',0,'d',0, 'c', 1/mu0, 'a', a_coeff, 'f', 0, 'Face', 2);

% generate mesh
generateMesh(model, 'Hmax', 0.002);

% solve PDE
result = solvepde(model);

The result of this solution is zero everywhere (not the expected solution). In addition, if I solve this in cartesian coordinates with a constant anonymous function for the coefficient (switching the definition of a_coeff), the simulation also returns zero everywhere. If I define a as zero in my call to specifyCoefficients, then the solution is correct.

Am I defining my coefficient function correctly? Is there a problem with my geometry (e.g. overlapping regions, if that is possible)?

[James Maniscalco]

I realize, incidentally, that my math is slightly wrong for the coefficient in cylindrical coordinates, but there should still be a non-zero solution (as evidenced by finding one in the PDE Modeler GUI).