The general result you are looking for (under the stated assumptions) looks like this:
For linear regression with $p$ predictor variables (you have two, $X$ and $X^2$) and an intercept, then with $n$ observations, $\mathbf{X}$ the $n \times (p+1)$ design matrix, $\hat{\beta}$ the $p+1$ dimensional estimator and $a \in \mathbb{R}^{p+1}$
$$ \frac{a^T\hat{\beta} - a^T \beta}{\hat{\sigma}
\sqrt{a^T(\mathbf{X}^T\mathbf{X})^{-1}a}} \sim t_{n-p-1}.$$
The consequence is that you can construct confidence intervals for any linear combination of the $\beta$ vector using the same $t$-distribution you use to construct a confidence interval for one of the coordinates.
In your case, $p = 2$ and $a^T = (0, x_2 - x_1, x_2^2 - x_1^2)$. The denominator in the formula above is the square root of what you compute as the estimate of the standard error (provided that this is what the software computes ...). Note that the variance estimator, $\hat{\sigma}^2$, is supposed to be the (usual) unbiased estimator, where you divide by the degrees of freedom, $n-p-1$, and not the number of observations $n$.
If I am understanding the question properly, you are asking about the non-linear least squares (NLS) model. Sticking to your notation, the NLS model is:
\begin{align}
y_n &= f(x_n,P) + \epsilon_n
\end{align}
Various regularity conditions are required on $f$, $X$, and $\epsilon$. The parameters are $P$, a $k$-vector. They are estimated by solving:
\begin{align}
min_P \sum_{n=1}^N \left(y_n-f(x_n,P) \right)^2
\end{align}
Assuming that you have found the global minimum and that it is interior to whatever the feasible set is for $P$ ($\mathbb{R}^k$, I guess?), the necessary FOC are:
\begin{align}
J_f\left(Y-f(X,P)\right)=0
\end{align}
The estimator, $P'$, defined as the interior solution to the minimization problem above and therefore solving the first-order condition immediately above, is consistent and asymptotically normal. It has an asymptotic variance of $\sigma^2_rH^{-1}$. The variance of the error in the original model is assumed to be $Var(\epsilon_n)=\sigma^2_r$, and it may be consistently estimated as $\hat{\sigma}_r^2=\frac{1}{N-k}\sum_{n=1}^N\left(y_n-f(x_n,P')\right)^2$.
These are standard results for the non-linear least squares model. My favorite reference for them is Amemiya, T (1985) Advanced Econometrics, section 4.3.
The upshot of all this is that, if you want a 95% confidence interval for, say, the $3^{rd}$ element of $P$, you can use:
\begin{align}
P_{3}' \pm 1.96\sqrt{\hat{\sigma}_r^2\left(H^{-1}\right)_{33}}
\end{align}
Notice that one of the regularity conditions is that all the variances of the $\epsilon$ are equal and that all of the covariances among elements of the $\epsilon$ are zero. The model in the link you provided contemplates that different $\epsilon_n$ have different variances (i.e. heteroskedasticity).
If that characterizes your application, then you need to modify your variance matrix. You can use a so-called sandwich estimator, like this:
\begin{align}
C &= \left( J_f^T J_f \right)^{-1}J_f^T \hat{\Sigma} J_f \left( J_f^T J_f \right)^{-1}\\
\hat{\Sigma} &= diag \left(\left(y_n-f(x_n,P') \right)^2 \right)
\end{align}
Then, the 95% confidence interval for the $3^{rd}$ element of $P$, would be:
\begin{align}
P_{3}' \pm 1.96\sqrt{C_{33}}
\end{align}
Further modifications would be required if there were to be correlations among the various elements of $\epsilon$.
Best Answer
No, it is not valid:
It can produce nonsense: for example if $a_1=-1$, $b_1=2$, $a_2=-2$ and $b_2=3$ then your statement would produce a confidence interval for the percentage change starting at -200% and increasing to -400%.
The width of the confidence interval ought to depend on any dependence between the two random variables.
Even if the the random variables are independent and positive, the confidence interval will depend on their distribution.