PDE Coefficients in the pde toolbox

Hi,

I'm going to solve the following elliptic equation in a domain which consists of two materials:

According to pde toolbox syntax the coefficients of the above (nonlinear) elliptic equation are:

c = K a = 0 f = d(K)/dz = d(K)/dh * dh/dz

but when I plot the flux terms (inside square brackets terms) using:

 uintrp = pdeintrp(p,t,u);
 [ux,uz] = pdegrad(p,t,u) ;
 for j=1:2
     idx = Sidx==j ;  % 'Sidx' contains the material index (1 or 2) for each triangle 
     qx(idx) = -F(j).K(uintrp(idx)) .* ux(idx) ;
     qz(idx) = -F(j).K(uintrp(idx)) .* (uz(idx)+1) ;
 end

 figure ; pdeplot(p,e,t,'flowdata',[qx; qz])

which should be conserved in the entire domain, I get this:

By setting "f=0" which is equivalent for:

the result is completely acceptable (the fluxes are equal everywhere). So it seems something is wrong with "f". For the reference "f" is calculated by:

 function f=fCoef(p,t,u,time)

 numElems = size(t,2) ;
 f = zeros(1,numElems) ;

 uintrp = pdeintrp(p,t,u); % Interpolated values at centroids
 [~,uz] = pdegrad(p,t,u); % Approximate derivatives

 for j=1:numElems    
    z1 = p(2,t(1,j)) ;
    z2 = p(2,t(2,j)) ;
    z3 = p(2,t(3,j)) ;
    zm = (z1+z2+z3)/3 ;

    if zm<=4
        f(j) = F(2).dKdh(uintrp(j)).* uz(j) ;
    else
        f(j) = F(1).dKdh(uintrp(j)).* uz(j) ;
    end
 end

Where am I missing something?

Thanks