When using bvp4c is it possible to specify von Neumann BCs for a first order equation

Hi,

I have a large set of ODEs to solve, but the problem lies in a temperature distribution. From a physical point I need a no flux condition, but the definition for the temperature gradient is

    dfdz(21)=-Tb*sum(dCb./Cb)-(delta*H)/(Cpmb*ub)*(Te-Tb)-delta/(ub*Cpmb)*rb'*dHr;

Is it possible to get bvp4c to use the derivative at the end directly, or do I need to do something ugly and inefficient like within the boundary condition function calculating

    -Tb*sum(dCb./Cb)-(delta*H)/(Cpmb*ub)*(Te-Tb)-delta/(ub*Cpmb)*rb'*dHr

And using that as my boundary condition?

If this is in fact the case, is there any way to pass parameters that depend on position to my boundary condition function, or would I need to recalculate these?

Many thanks, Doug