NumericQ prevents evaluation at a numerical value sometimes

Question

I set up a symbolic integral, to be evaluated only when some of the symbols have been replaced by numbers. The evaluation after assigning numbers to the symbols is still symbolic for some reason. The code is (after clearing all)

$Assumptions=0 < r < 1 && 0 < ah < al && 0 < s && t \[Element] Reals && llow < lhi
rev[t_] := Exp[t]/(1 + Exp[t])
yh1[t_] := \[ExponentialE]^(1/2 t (-1 - Sqrt[1 + 8 r s^2]))
yh2[t_] := \[ExponentialE]^(1/2 t (-1 + Sqrt[1 + 8 r s^2]))
(* Wronskians *)
wrh[t_] := yh1[t]*yh2'[t] - yh1'[t]*yh2[t]
(* Particular solutions of the inhomogeneous equations *)
uh[t_] := rev[t] - ah
inth1[t_, llow_?NumericQ, r_?NumericQ, s_?NumericQ, ah_?NumericQ] := 
 Integrate[yh2[z]*uh[z]/wrh[z], {z, llow, t}]
inth2[t_, llow_?NumericQ, r_?NumericQ, s_?NumericQ, ah_?NumericQ] := 
 Integrate[yh1[z]*uh[z]/wrh[z], {z, llow, t}]
yhp[t_, llow_, r_, s_, ah_] := -yh1[t]*inth1[t, llow, r, s, ah] + 
 yh2[t]*inth2[t, llow, r, s, ah] 
yhp[t, llow, r, s, ah]
% /. {llow -> -2, lhi -> 2, ah -> 0.1, r -> 1, s -> 1}

The output of yhp[t, llow, r, s, ah] is

-\[ExponentialE]^(1/2 t (-1 - Sqrt[1 + 8 r s^2])) inth1[t, llow, r, s,ah] + 
 \[ExponentialE]^(1/2 t (-1 + Sqrt[1 + 8 r s^2]))   inth2[t, llow, r, s, ah] 

like it's supposed to be. But the output of % /. {llow -> -2, lhi -> 2, ah -> 0.1, r -> 1, s -> 1} has r,s etc still in it.

I tried to reproduce the problem in a simple example, but the example works fine. It is

Clear[x, k, int, f]
$Assumptions = {x, k} \[Element] Reals
int[x_, k_?NumericQ] := 
 Integrate[Exp[z]/((1 + Exp[z])*Exp[k*z]), {z, 0, x}]
f[x_, k_] := x*int[x, k]
f[x, k]
% /. {k -> 2}

and the final output is

x (1 - \[ExponentialE]^-x - Log[2] + Log[1 + \[ExponentialE]^-x]) 

as it's supposed to be.

Answer 1

Your functions depend on many variables that have not been included in their definitions. This is a reformulation that seems to work for me:

$Assumptions = 0 < r < 1 && 0 < ah < al && 0 < s && 
     Element[t, Reals] && llow < lhi;
rev[t_] := Exp[t]/(1 + Exp[t])
yh1[t_, r_, s_] := E^((1/2)*t*(-1 - Sqrt[1 + 8*r*s^2]))
yh2[t_, r_, s_] := E^((1/2)*t*(-1 + Sqrt[1 + 8*r*s^2]))
wrh[t_, r_, s_] := yh1[t, r, s]*Derivative[1, 0, 0][yh2][t, r, 
         s] - Derivative[1, 0, 0][yh1][t, r, s]*yh2[t, r, s]
uh[t_, ah_] := rev[t] - ah
inth1[t_, (llow_)?NumericQ, (r_)?NumericQ, (s_)?NumericQ, 
     (ah_)?NumericQ] := NIntegrate[(yh2[z, r, s]*uh[z, ah])/
       wrh[z, r, s], {z, llow, t}]
inth2[t_, (llow_)?NumericQ, (r_)?NumericQ, (s_)?NumericQ, 
     (ah_)?NumericQ] := NIntegrate[(yh1[z, r, s]*uh[z, ah])/
       wrh[z, r, s], {z, llow, t}]
yhp[t_, llow_, r_, s_, ah_] := 
   (-yh1[t, r, s])*inth1[t, llow, r, s, ah] + 
     yh2[t, r, s]*inth2[t, llow, r, s, ah]
yhp[t, llow, r, s, ah] /. {t -> 2, llow -> -2, lhi -> 2, 
     ah -> 0.1, r -> 1, s -> 1}
(* 3.02156 *)