Can anyone tell me how I can use `NDSolveValue` to model reflective and constant heat flux boundary conditions. I am solving the heat equation. Essentially I have a PDE that is dependent on time, radius, and axial length. I want to solve it such that the spatial derivative of the temperature is equal to 0 on one boundary and equal to a constant on another boundary. Can anyone explain how I can set up `NDSolveValue` with these types of boundary conditions? I don't think it can be achieved using [`NeumannValue`][1] because that just sets the entire differential equation to a value. I have tried a lot of different approaches but nothing seems to be working for me. Can anyone please recommend a way to achieve this?

Edit: Here is some code

    tf = 50; Ti = 100; 
    Ls =  250; Lito = 5;  Lsl = 230;
    Ltot = Ls + Lito + Lsl;
    R = 1500; k = 1; 
    
    eqn = r*\!\(
    \*SubscriptBox[\(\[PartialD]\), \(t\)]\(T[t, r, z]\)\) - r*k*\!\(
    \*SubscriptBox[\(\[PartialD]\), \(z, z\)]\(T[t, r, z]\)\) - k*\!\(
    \*SubscriptBox[\(\[PartialD]\), \(r\)]\(T[t, r, z]\)\) - k*r*\!\(
    \*SubscriptBox[\(\[PartialD]\), \(r, r\)]\(T[t, r, z]\)\);
    Subscript[\[CapitalGamma], 
      D] = {DirichletCondition[T[t, r, z] == Ti, z == Ltot], 
       DirichletCondition[T[t, r, z] == Ti, r == R]};
    
    BCr = NDSolveValue[{eqn == 0, Subscript[\[CapitalGamma], D], \!\(
    \*SubscriptBox[\(\[PartialD]\), \(r\)]\(T[t, 0, z]\)\) == 0, \!\(
    \*SubscriptBox[\(\[PartialD]\), \(z\)]\(T[t, r, 0]\)\) == -100, 
        T[0, r, z] == Ti}, T, {t, 0, tf}, {r, 0, R}, {z, 0, Ltot}];

  [1]: http://reference.wolfram.com/language/ref/NeumannValue.html