problem in calculating flux in pde toolbox

Hi everybody,

I'm trying to calculate and plot the flux terms in the following equation in the pde toolbox:

fluxes are the inside square bracket terms. According to pde toolbox syntax, pde coefficients become:

c = K

a = 0

f = d(K)/dz = d(K)/dh * dh/dz

when I try to plot the fluxes by these lines the result is incorrect (violence of mass conservation law):

 uintrp = pdeintrp(p,t,u) ;
 [ux,uz] = pdegrad(p,t,u) ;
 qx = -K(uintrp) .* ux ;
 qz = -K(uintrp) .* (uz + 1) ;
 pdeplot(p,e,t,'flowdata',[qx; qz])

But surprisingly when I change the 'qz' as follows everything seems OK!:

 qz = -K(uintrp) .* (uz - 1) ;

so I'm a bit confused since the flux in the main equation has '+' in it not '-'. What am I missing here?

thanks,

[Mohammad Monfared]

The problem was simply because of a '-' in a function evaluating dK/dh !