What if the resulted error is so large? Suppose that we need to do the error analysis for the quantity:
$$F = \frac{1}{A-B}$$
with $A=1.08\pm0.02$ and $B = 1.05\pm0.03$. As usual, it could be evaluated by
$$\bar{F} = \overline{A}-\overline{B} \simeq 33.333$$
$$\delta F = \bar{F}\frac{\delta (A-B) }{\overline{A-B}} = \bar{F}\times\frac{\sqrt{(\delta A)^2 +(\delta B)^2} }{\overline{A}-\overline{B}} \simeq 40.062$$
So, F = $33.33 \pm 40.06$, which looks ridiculous. So, is this error analysis problematic?
(In the above example, A is close to B but always larger than B, at least theoretically)
 A: This is a case in which the standard statistical tool of Gaussian error propagation fails.
When you write a measured quantity as "$\text{quantity} = \text{value} \pm \text{error}$" this is typically interpreted as modelling it with a normal distribution, $\text{quantity} \sim \mathcal{N}(\text{quantity}; \mu = \text{value}, \sigma = \text{error})$.
The formula you quote for error propagation is derived assuming this, as well as considering only the linear order in the change of the quantity.
If the errors were much smaller this would work decently well: say, if you had $A = 1.08 \pm 10^{-6}$ and $B = 1.06 \pm 10^{-6}$, you'd get
a mean of $1/(A-B) = 50$ and an error of $50^2 \sqrt{2} 10^{-6} \approx 0.0035$, both of which agree with what we get if we simulate two normally-distributed populations: the python code is
import numpy as np
N = 10_000
rng = np.random.default_rng(seed=1)
A = rng.normal(loc=1.08, scale=1e-6, size=N)
B = np.random.normal(loc=1.06, scale=1e-6, size=N)
print(f'mean = {np.mean(1/(A-B)):.0f}')
print(f'error = {np.std(1/(A-B)):.4f}')

mean = 50
error = 0.0035

This works here because the function $1/(A-B)$ in the region of interest, basically $A-B \sim 0.02 \pm \text{a few times} \sqrt{2} \times 10^{-6}$, is almost perfectly linear. In this limit, the formula you quoted works well.
Why does this not extend to the values you quoted for the errors?
The region of interest now includes a place where $A-B$ goes to zero and becomes negative: $1/(A-B)$ is very much not a linear function here, and cannot be approximated as such.
You're getting errors on the order of a few times 10; if you try the simulation approach you'll get something like
import numpy as np
N = 10_000
rng = np.random.default_rng(seed=1)
A = rng.normal(loc=1.08, scale=.02, size=N)
B = rng.normal(loc=1.06, scale=.03, size=N)
print(f'mean = {np.mean(1/(A-B)):.0f}')
print(f'error = {np.std(1/(A-B)):.4f}')

mean = 27
error = 2169.0462


which is completely different!
This is not even the full picture: the actual number for the error is meaningless, in that it very much depends on the statistics of the simulation; increasing $N$ to $10^{5}$, for example, makes it increase to 25183.
This is an indication that we're not modeling the problem correctly: indeed, if you try to do the calculation analytically you'll find, as @Farcher noted in a comment, that the mean and variance of this distribution do not exist: the integrals which define them diverge.
So, the problem as stated is ill-posed: basically, you cannot really say much about the distribution of $1/(A-B)$, except for the fact that it can take on very high values.
If this is indeed the true problem statement, perhaps the best you can do is give a constraint on the median of the distribution.
You mentioned the fact that you have the constraint $A>B$; this does not really fix the problem, since you still get the divergence when $A-B$ approaches zero.
On the other hand, if you have some physical reason based on your measurement apparatus to believe that the difference $A-B$ cannot really approach 0 but must instead be larger than some value, you'll get a much more interesting result - which, however, will heavily depend on that lower bound, since it's there that the quantity $F$ can become really large.
Still, you'll need to account for the nonlinearity of the problem, and the formula you quoted will not do.
A: 
$F = \frac{1}{A-B}$


with $A=1.08\pm0.02$ and $B = 1.05\pm0.03$.


...which looks ridiculous. So, is this error analysis problematic?

It's not ridiculous. It might be problematic depending on how well you need to determine $F$.
But, given the uncertainties you provided, your $F$ could turn out to be infinite (given the error bars). E.g., $\frac{1}{1.08 - 1.08}$. So, it is no wonder the resulting errors are large.

(In the above example, A is close to B but always larger than B, at least theoretically)

Unfortunately, this is not the case for your measurements, given your error bars. Within the error we could have $A=1.06$ and $B=1.08$.
A: I think that this is not be unreasonable. F is a large number generated from 2 small numbers on a hyperbolic curve. As you get close to the y axis, the curve goes essentially vertical.
According to Taylor, "Introduction to Error Analysis": $$ \delta F = \sqrt{\left(\frac{\partial F}{\partial A}\delta A\right)^{2}+\left(\frac{\partial F}{\partial B}\delta B\right)^{2}}$$
This formula is essentially the same as your calculation and gives the same results. The difference is this form reveals the dependency on the rate of change of F with A and with B.
$$\frac {\partial F}{\partial A} = \frac {-1}{\left(A-B\right)^2} = \frac {-1}{F^2} = 1111$$
It may seem ridiculous at first look, but look again.
The other side of the coin is trying to generate a small number from two large numbers that appear to have reasonable uncertainty but the small quantity is swamped by noise.
