I have been trying to prove that origin $(\theta,k) = 0$ of the non-autonomous system $$ \begin{aligned} \dot{\theta}(t) &= -b \theta(t) + k(t) u \\ \dot{k}(t) &= \epsilon \text{ } h(\theta(t),k(t),u)\end{aligned}$$ with sufficiently small $1 \gg \epsilon > 0$, $b > 0$, Lipschitz continuous $h : \mathbb{R}^3 \to \mathbb{R}$ and $u \neq 0$ is locally stable, if $h$ satisfies $$\rm{sign}[h(\theta,k,u)] = -\rm{sign}[\theta \text{ } u].$$
I can't show that this is true for such a seemingly simple problem. From numerical simulations, it is clearly true regardless of h. However, showing it theoretically has been a pain. Can't apply linearization because have no information about the derivative of $h$ with respect states. Can't find a Lyapunov function that will reveal the stability since there is no explicit form of $h$.
What I think I can show is : if the $h$ condition is not satisfied in any open set containing the origin, then the origin cannot be stable. However, this doesn't show the claim.
Is there a way to show the claim?
Because of your small gain $\epsilon$ one can use time scale separation, such that from the perspective of the dynamics of $k(t)$ the variable $\theta(t)$ seems to follow a quasi-steady-state of
$$ \bar{\theta}(t) = \frac{k(t)\,u}{b}. \tag{1} $$
The error between the actual $\theta(t)$ and $(1)$ is defined as $\tilde{\theta}(t) = \theta(t) - \bar{\theta}(t)$ and can be shown to have to following dynamics
$$ \dot{\tilde{\theta}}(t) = -b\,\tilde{\theta}(t) - \frac{\epsilon\,h(\theta,k,u)\,u}{b}. \tag{2} $$
Using the convolution integral it can be shown that from $(2)$ it follows that
$$ \tilde{\theta}(t) = e^{-b\,t}\,\tilde{\theta}(0) + \epsilon \int_0^t e^{b\,(\tau-t)}\,\frac{h(\theta(\tau), k(\tau), u)\,u}{b}\,d\tau. \tag{3} $$
I am not too familiar with giving complete proofs when it comes to time scale separation, so this might not be completely correct but hopefully at least gives some intuition and insights in how to proceed. Namely, $(2)$ should be BIBO, so if $h(\theta,k,u)$ remains bounded then that integral should also be bounded. Therefore, after enough time $t_0$, in which the transient $e^{-b\,t}\,\tilde{\theta}(0)$ has practically vanished, it should hold that
$$ \theta(t) = \bar{\theta}(t) + O(\epsilon), \tag{4} $$
where $O(\epsilon)$ can be seen as a perturbation. For the unperturbed system one can substitute $(1)$ for $\theta(t)$ in the sign constraint, which yields
$$ \text{sign}[h(\bar{\theta}(t), k(t), u)] = -\text{sign}[\bar{\theta}(t)\,u] = -\text{sign}\!\left[\frac{k(t)\,u^2}{b}\right], \tag{5} $$
and since $u \neq 0$ and $b > 0$ it follows that $(5)$ is equivalent to
$$ \text{sign}[h(\bar{\theta}(t), k(t), u)] = -\text{sign}[k(t)]. \tag{6} $$
Hopefully it is clear that from $(6)$ for the unperturbed system it follows that $k(t)$ and thus $\theta(t)$ would be asymptotically stable, since $\text{sign}[\dot{k}(t)] = -\text{sign}[k(t)]$.
The perturbed sign constraint would give
$$ \text{sign}[h(\theta(t), k(t), u)] = -\text{sign}\!\left[\frac{k(t)\,u^2}{b} + u\,O(\epsilon)\right], \tag{7} $$
so once $k(t)$ becomes small enough the contribution of $u\,O(\epsilon)$ might become significant. Therefore, I am not entirely sure how to proceed here.
I also tried to look for a potential Lyapunov function for $t > t_0$. Using $\tilde{\theta}(t)$ and $(2)$ a potential candidate might be
$$ V(\tilde{\theta}(t), k(t)) = \frac{1}{2} \epsilon\,\tilde{\theta}(t)^2 + \frac{1}{2} k(t)^2, \tag{8} $$
for which it can be shown that it has the following time derivative
$$ \dot{V}(\theta(t), k(t)) = \epsilon \left( -b\,\tilde{\theta}(t)^2 - \tilde{\theta}(t) \frac{\epsilon\,h(\theta,k,u)\,u}{b} + k(t)\,h(\theta,k,u) \right). \tag{9} $$
The term $\tilde{\theta}(t)\,\epsilon\,h(\theta,k,u)\,u / b$ should be negligibly, since it is multiplied by $\epsilon$, but this still leaves the issue that $k(t)\,h(\theta,k,u)$ might not be negative definite in $k(t)$, even with the assumption that $t>t_0$.
Another approach could be to use the circle criterion, since a sign constraint would translate nicely to the sector bound $[0,\infty]$. However, in order for this approach to work it is required to have the following tighter sector constraint
$$ h(\theta,k,u)\,[\epsilon\,h(\theta,k,u) - K\,\theta\,u] \leq 0, $$
with $0 < K < \frac{b^2}{u^2}$. So if $\epsilon$ can chosen to be arbitrarily close to zero then this constraint essentially requires that $h(\theta,k,u)$ is continuous at $\theta=0$. If you would be able to impose this constraint and would like more details regarding to this approach feel free to ask.