In short, is there a clean formula for $\frac{d^n}{dx^n}\operatorname{sech}^2(x)?$
For convenience, let
\begin{align} f(x)&=\operatorname{sech}^2(x),\\\\ g(x)&=\tanh(x). \end{align}
Then
\begin{align} f'(x)&=-2\operatorname{sech}^2(x)\tanh(x)& &=af(x)g(x),\\\\ g'(x)&=\operatorname{sech}^2(x)& &=f(x), \end{align}
where $a=-2.$ We also have $g^2(x)=1-f(x).$
Hence $f^{(n)}(x)$ can always be written as a sum where each term is of the form $c_{k,p} [f(x)]^k[g(x)]^p$ where $p=0$ or $1$. So the question becomes: is there a formula for the $c_{k,p}$'s?
I would expect some pattern to emerge from repeated uses of the power rule combined with the recursion of the above formulas, but I'm not sure how to proceed. Also, I know sometimes coefficients of special polynomials can encode information about $n^\text{th}$ derivatives (e.g., Hermite polynomials and $e^{-x^2}$), and I would find it acceptable to express the answer in terms of such coefficients if possible.
Using the General Leibniz rule, $$\mathrm D^n\operatorname{sech}^2=\sum_{k=0}^n\binom{n}{k}~\mathrm D^{n-k}\operatorname{sech}~\cdot~\mathrm D^{k}\operatorname{sech}$$ So our problem reduces to finding $\mathrm D^n \operatorname{sech}$.
To clean up the notation a bit, let $f_n(\cdot,\cdot)$ be a polynomial such that $$\mathrm D^n(\operatorname{sech})(t)=f_n(\operatorname{sech}(t),\tanh(t))$$ The first few of them are $$ \begin{array}{l} f_{0}( x,y) =x\\ f_{1}( x,y) =-xy\\ f_{2}( x,y) =-x^{3} +xy^{2}\\ f_{3}( x,y) =5x^{3} y-xy^{3}\\ f_{4}( x,y) =5x^{5} -18x^{3} y^{2} +xy^{4}\\ f_{5}( x,y) =-61x^{5} y+58x^{3} y^{3} -xy^{5}\\ f_{6}( x,y) =-61x^{7} +479x^{5} y^{2} -179x^{3} y^{4} +xy^{6}\\ f_{7}( x,y) =1385x^{7} y-3111x^{5} y^{3} +543x^{3} y^{5} -xy^{7}\\ f_{8}( x,y) =1385x^{9} -19\ 028x^{7} y^{2} +18\ 270x^{5} y^{4} -1636x^{3} y^{6} +xy^{8} \end{array}$$ Hmm. Some patterns are appearing in the first and last terms, but the middle terms are a bit murky. Still, it is clear that only $\operatorname{sech}$ and $\tanh$ will be present in the result. With that in mind it is pertinent to try to evaluate $$\mathrm D(\operatorname{sech}^l\tanh^m)$$ Where the superscript denotes exponentiation, not composition. One finds $$\mathrm D(\operatorname{sech}^l\tanh^m)=m\operatorname{sech}^{2+l}\tanh^{m-1}-l\operatorname{sech}^{l}\tanh^{m+1}$$ So in our polynomials, terms of the form $x^ly^m$ in $f_n(x,y)$ change into $mx^{2+l}y^{m-1}-lx^ly^{m+1}$ in $f_{n+1}(x,y)$. IMPORTANT: This also makes it clear that $f_n(x,y)$ will only contain terms of the form $x^{n+1-k}y^k$. This suggests the amazing recurrence relation $$f_{n}(x,y)=x^2~\partial_y f_{n-1}(x,y)-xy~\partial_xf_{n-1}(x,y)$$ So, letting $$f_n(x,y)=\sum_{k=0}^{n+1}a_{n,k}x^{n+1-k}y^k$$ We have $$\sum_{k=0}^{n+1}a_{n,k}x^{n+1-k}y^k=x^2\sum_{k=0}^{n}a_{n-1,k}~k~x^{n-k}y^{k-1}-xy\sum_{k=0}^{n}a_{n-1,k}(n-k)x^{n-k-1}y^k$$ In the first sum we now change index $k= j+1$ and in the second sum we change $k=l-1$ to obtain $$\sum _{k=0}^{n+1} a_{n,k} x^{n+1-k} y^{k} $$ $$=x^{2}\sum _{j=-1}^{n-1} a_{n-1,j+1}( j+1) x^{n-j-1} y^{j} -xy\sum _{l=1}^{n+1} a_{n-1,l-1}( n-l+1) x^{n-l} y^{l-1}$$ Removing the first summand from the first sum (since it's zero), renaming our indices, and distributing the outer factors we get
$$\sum _{k=0}^{n+1} a_{n,k} x^{n+1-k} y^{k} =\begin{matrix} \sum _{k=0}^{n-1} a_{n-1,k+1}( k+1) x^{n+1-k} y^{k}\\ -\sum _{k=1}^{n+1} a_{n-1,k-1}( n-k+1) x^{n+1-k} y^{k} \end{matrix}$$ Adding and subtracting some terms we get $$\sum _{k=0}^{n+1} a_{n,k} x^{n+1-k} y^{k} =\begin{matrix} -a_{n-1,n+1}( n+1) xy^{n} -a_{n-1,n+2}( n+2) y^{n+1}\\ +\sum _{k=0}^{n+1} a_{n-1,k+1}( k+1) x^{n+1-k} y^{k}\\ +a_{n-1,-1}( n+1) x^{n+1}\\ -\sum _{k=0}^{n+1} a_{n-1,k-1}( n-k+1) x^{n+1-k} y^{k} \end{matrix}$$ But for obvious reasons, we must have $$a_{n-1,n+1}=a_{n-1,n+2}=a_{n-1,-1}=0$$ So in all we have $$\sum _{k=0}^{n+1} a_{n,k} x^{n+1-k} y^{k} =\sum _{k=0}^{n+1}( a_{n-1,k+1}( k+1) -a_{n-1,k-1}( n-k+1)) x^{n+1-k} y^{k}$$ And so by a linear algebra argument we have a recurrence relation $$\boxed{a_{n,k}=(k+1)a_{n-1,k+ 1}-(n-k+1)a_{n-1,k-1}}$$ With an initial conditions $a_{0,0}=1,a_{0,1}=0$ and boundary conditions $$a_{n,k}=0\text{ if }k<0\text{ or }k>n+1$$
A quick test: we know that for example we should have $$a_{1,0}=0,a_{1,1}=-1,a_{1,2}=0$$ We have $$a_{1,0}=a_{0,1}-(1-0+1)a_{0,-1}=0$$ $$a_{1,1}=2a_{0,2}-(1-1+1)a_{0,0}=-1$$ $$a_{1,2}=3a_{0,3}-(1-2+1)a_{0,1}=0$$ Our formula is working as expected. So using this definition of $a_{n,k}$ we write $$\mathrm{D}^n(\operatorname{sech})=\sum_{k=0}^{n+1}a_{n,k}\operatorname{sech}^{n+1-k}\tanh^k$$ Therefore finally,
$$\mathrm{D}^n(\operatorname{sech}^2)=\sum_{k=0}^n\left[\binom{n}{k}\left(\sum_{j=0}^{n-k+1}a_{n-k,j}\operatorname{sech}^{n-k+1-j}\tanh^{j}\right)\left(\sum_{l=0}^{k+1}a_{k,l}\operatorname{sech}^{k+1-l}\tanh^l\right)\right]$$ I don't think we can simplify this any further!