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.
