Maximize the variance of a function of random variable

I have a function $f(X)=\exp\left(\frac{-\gamma^2}{a^2X+b^2}\right)$ where $X \sim \mathrm{Binomial}(n,p)$. I am interested in finding the value of $\gamma^2$ which maximizes the variance of $f(X)$ for given values of $a$ and $b$.

Definition of variance: $$ \operatorname{var}(f(X))=\mathrm{E}\left(f^2(X)\right)-\left(\mathrm{E}\bigl(f(X) \bigl) \right)^{2} $$ From Expectation of Function of Discrete Random Variable:

$$E\left(f^2(X)\right)=\sum_{k = 0}^{n} \exp\left(\frac{-2\gamma^2}{a^2k+b^2}\right) {n \choose k} p^{k} (1-p)^{n-k}$$ $$E\left(f(X)\right)=\sum_{k = 0}^{n} \exp\left(\frac{-\gamma^2}{a^2k+b^2}\right) {n \choose k} p^{k} (1-p)^{n-k}$$

To me, it seems NOT straightforward how to solve $\mathrm{d}[\operatorname{var}(f(X))]/\mathrm{d}\gamma=0$ to find an expression for optimal $\gamma^2$.

Hence, I am thinking of at least finding a numerically accurate solution through curve-fitting.

My aim is to come up with a "working formula" for optimal $\gamma$ in terms of the other parameters.

Also, I have MATLAB and I am wondering how formulate this problem and find an expression using MATLAB as well.

Some context on $f(X)$: $f(X)$ comes from the complementary CDF of an exponential distribution. There are $X$ signal components with power $a^2$ and a non-signal component with power $b^2$. Also, $\gamma^2$ is a threshold at which we evaluate the function.

This answer tries to come up with a working formula, although an approximate one.

I use the following approximate value to the variance of a function of a random variable: $$\operatorname{Var}\left[f(X)\right] \approx \left(f'(\operatorname{E}\left[X\right])\right)^2\operatorname{Var}\left[X\right] $$

In this case, $\operatorname{E}\left[X\right]=np$ and $$\frac{d}{d X}f|_{X=np}=\exp \left(\frac{-\gamma^2}{a^2np+b^2}\right)\frac{a^2\gamma^2}{(a^2np+b^2)^2}$$

I set aside the constant $\operatorname{Var}\left[X\right]$.

Maximizing the square of a function is equivalent to maximizing the function itself.

Letting $$ z:= \frac{\gamma^2}{a^2 n p+b^2}$$ $$ C:= \frac{a^2}{a^2np+b^2}$$ the next step is to maximize the function $z \mapsto C z \exp(-z)$ over the positive reals.

The standard procedure gives $z=1$ which amounts to $$\gamma^2= a^2 np+b^2$$.

I hope this helps.