1. Can we numerically detect stiffness just by applying explicit methods?
Suppose you have an initial value problem for some ODE on $[0,10]$.
You take considerably large stepsize $\tau=1$ and an explicit Euler
method, make your calculations with constant step size $\tau$ and
get these points:
You estimate the error and it appears to be big. Okay, then you take
$\tau=0.1$ and obtain
The error estimate is acceptable now. Stepsize $\tau=0.1$ is small
relatively to $[0,10]^\star$.
So, is the problem stiff? The answer is NO! Small stepsize here is required to correctly reproduce the oscillations of solution.
The problem we've solved is $$ y'(t)=-2\cos \pi t,\quad y(0)=1. $$
Now you take another ODE on the same interval, $\tau=1$ and explicit Euler gives you almost the same numerical solution:
You take $\tau=0.1$ and now the numerical solution is
The error estimate now is small. Stepsize $\tau=0.1$ is small relatively to $[0,10]^\star$.
Is this problem stiff? YES! We've made very much small steps to reproduce the solution which is changing very slowly. This is irrational! The magnitude of time step here is limited by the stability properties of explicit Euler.
This problem is
y'(t)=-2 y(t)+\sin t/2,\quad y(0)=1.
$^\star$ Note that the number of step sizes can be much greater if we take longer interval of integration.
Conclusion: the information about timesteps and corresponding errors is not sufficient to detect stiffness. You should also look at the obtained solution. If it varies slowly and stepsize is very small, the problem is most likely to be stiff. If the solution oscillates rapidly and you trust your error estimation technique then this problem is not stiff.
2. How to determine the maximum stepsize which allows to integrate stiff problem with explicit method?
If you use some black-box explicit solver with automatic step control then you need nothing to do: the software will take the required stepsize adaptively.
But suppose you want to get the solution manually with constant stepsize, or just want to estimate how many hours you should wait until your explicit method crunches the problem.
Then you should know the spectrum of your Jacobi matrix. Suppose it is real and lies in $[\Lambda,0]$ (in your example $\Lambda=-1000$).
Then you should calculate the real stability interval (stability domain in complex case) of your explicit method. It is not too difficult, you need just consult any textbook on this topic. In the case of explicit Euler this interval is $[-2,0]$. Now, if you want your solution not to blow-up, you should take $\tau$ such that $\Lambda\tau$ lies in the stability interval, i. e. in our case
If you want more consistency, you should take
since for $1/|\Lambda|<\tau \leq 2/|\Lambda|$ your solution will most likely produce unnatural (but fading) oscillations.
Of course such an analysis is mostly applicable for linear problems with known spectrum. For more practical problems we should rely on numerical methods of stiffness detection (see references and comments in other answers).