The variance of the sum of two variables must be calculated with a term accounting for the covariance of those two variables.
$$
Var(aX + bY) = a^2Var(X) + b^2Var(Y) + 2ab Cov(X,Y)
$$
Note that the coefficients on the variables are also squared in the first two terms of that equation. If X and Y are fully independent then the third term is zero of course.
Recall the expression for the posterior variance of a Gaussian Process: $\operatorname{var}_{n}\left(f\left(\mathbf{x}_{*}\right)\right)=k\left(\mathbf{x}_{*}, \mathbf{x}_{*}\right)-\mathbf{k}_{*}^{\top}\left(K_n+\sigma^{2} I_n\right)^{-1} \mathbf{k}_{*}$. We want to show that $\text{var}_{n-1}(f(x_*)) \geq \text{var}_n(f(x_*))$. The first term is the same for both $\text{var}_n$ and $\text{var}_{n-1}$, so we just need to break down the second term in a way such that we can compare between $\text{var}_n$ and $\text{var}_{n-1}$.
Following the hint to look at section A.3 in the book, we see that for the partitioned matrix $A$ and its inverse $A^{-1}$,
\begin{align}
A = \begin{pmatrix} P & Q \\ R & S \end{pmatrix}, \qquad A^{-1} = \begin{pmatrix} \tilde P & \tilde Q \\ \tilde R & \tilde S \end{pmatrix},
\end{align}
\begin{align}
\tilde P & = P^{-1} + P^{-1}QMRP^{-1},\\
\tilde Q & = -P^{-1}QM,\\
\tilde R & = -MRP^{-1},\\
\tilde S & = M,\\
\text{where}\ M & = (S - RP^{-1}Q)^{-1}.\\
\end{align}
We can decompose the Gram matrix with added noise $K_n + \sigma^2I_n$ as follows:
\begin{align}
K_n + \sigma^2 I_n =
\begin{pmatrix}
K_{n-1} + \sigma^2 I_{n-1} & k_{n-1}(x')\\k_{n-1}(x')^\top & k(x',x') + \sigma^2
\end{pmatrix}
\end{align}
where $x'$ is the $n$th point sampled.
Writing out the decomposed inverse, we get
\begin{align}
(K_n + \sigma^2 I_n)^{-1} = \begin{pmatrix}\kappa + \kappa k_{n-1}(x')Mk_{n-1}(x')^\top\kappa & -\kappa k_{n-1}(x')M \\ -Mk_{n-1}(x')^\top \kappa & M \end{pmatrix},
\end{align}
with
\begin{align}
M & = (k(x', x') + \sigma^2 - k_{n-1}(x')^\top(K_{n-1} + \sigma^2I_{n-1})^{-1}k_{n-1}(x'))^{-1}, \\
\kappa & = (K_{n-1} + \sigma^2I_{n-1})^{-1}
\end{align}
Finally we can actually compute the term $k_*^\top(K_n + \sigma^2I_n)^{-1}k_*$ by multiplying through:
\begin{align}
k_*^\top(K_n + \sigma^2I_n)^{-1}k_* = & \begin{pmatrix}k_{n-1}(x^*)\\k'(x^*)\end{pmatrix}^\top\begin{pmatrix}\kappa + \kappa k_{n-1}(x')Mk_{n-1}(x')^\top\kappa & -\kappa k_{n-1}(x')M \\ -Mk_{n-1}(x')^\top\kappa & M \end{pmatrix}\begin{pmatrix}k_{n-1}(x^*)\\k'(x^*)\end{pmatrix}\\
& = \begin{pmatrix} k_{n-1}^\top(x^*)\kappa + k_{n-1}^\top(x^*)\kappa k_{n-1}(x')Mk_{n-1}(x')^\top\kappa - k'(x^*)Mk_{n-1}(x')^\top \kappa \\ -k_{n-1}^\top(x^*)\kappa k_{n-1}(x')M + k'(x^*)M \end{pmatrix} ^\top
\begin{pmatrix}k_{n-1}(x^*)\\k'(x^*)\end{pmatrix}
\end{align}
where $\kappa = (K_{n-1} + \sigma^2I_{n-1})^{-1}$ and $k'(x^*) = k(x', x^*)$.
We end up with
\begin{align}
k_*^\top(K_n + \sigma^2I_n)^{-1}k_* & = k_{n-1}^\top(x^*)\kappa k_{n-1}(x^*) + k_{n-1}^\top(x^*)\kappa k_{n-1}(x')Mk_{n-1}(x')^\top\kappa k_{n-1}(x^*)\\
& - k'(x^*)Mk_{n-1}(x')^\top\kappa k_{n-1}(x^*) - k_{n-1}(x^*)^\top\kappa k_{n-1}(x')Mk'(x^*) + k'(x^*)Mk'(x^*).
\end{align}
Note that the first term is simply the relevant term from $\text{var}_{n-1}$, so we just need to show that the subsequent terms are non-negative in order to show that the term is larger for $n$ compared to $n-1$ (and so the variance is smaller). Note that $M$ is the reciprocal of the variance at $x'$ plus $\sigma^2$, so it's positive and we can factor it out.
We are left with
\begin{align}
\alpha^2 - 2k'(x^*)\alpha + k'(x^*)^2,
\end{align}
where $\alpha = k_{n-1}^\top(x^*)\kappa k_{n-1}(x')$ (a scalar). So we have
\begin{align}
k_*^\top(K_n + \sigma^2I_n)^{-1}k_* & = k_{n-1}^\top(x^*)\kappa k_{n-1}(x^*) + (\alpha - k'(x^*))^2\\
& = k_{n-1}(x^*)^\top(K_{n-1} + \sigma^2I_{n-1})^{-1}k_{n-1}(x^*) + \tfrac{1}{M}(\alpha - k'(x^*))^2\\
& \geq k_{n-1}(x^*)^\top(K_{n-1} + \sigma^2I_{n-1})^{-1}k_{n-1}(x^*).
\end{align}
So the variance after $n$ points is smaller than the variance after $n-1$ points, with equality achieved when the quantity $\tfrac{1}{M}(k_{n-1}^\top(x^*)(K_{n-1} + \sigma^2I_{n-1})^{-1} k_{n-1}(x') - k(x^*, x'))$ is zero.
Best Answer
Can you not use $e^x \approx 1+ x+ \frac{X^2}{2} $ then $Var(e^X) \approx Var(X)+ \frac{Var(X^2)}{2} $ ?
It seems to me you want 2nd order approximation ..