The heat diffusion equation on some domain $\Omega$ with Neumann boundary conditions on $\partial\Omega$ and normal $n$ is given as: \begin{alignat}{3} \partial_t u(t,x) &= \Delta u(t,x), &\quad& x\in\Omega, \, t\in(0,\infty],\\ \partial_n u(t,x) &= 0, &\quad& x\in\partial\Omega, \, t\in(0,\infty], \\ u(0,x) &= f(x), &\quad& x\in\Omega. \end{alignat} A discretisation in space with $L\approx -\Delta$, $u \approx w:[0,\infty] \to \mathbb{R}^n$, and $f\approx b\in\mathbb{R}^n$ leads to a system of ODEs: $$d_t w(t) = -Lw(t), \quad w(0) = b.$$ The solution is then given as $$w(t) = \exp(-t L)b = V \begin{bmatrix} \exp(-t\lambda_1) && \\ &\ddots& \\ && \exp(-t\lambda_n) \end{bmatrix} V^T, \quad L = V\Lambda V^Tb.$$ This formulation has no issues whatsoever with $t\to\infty$ and leads to the average of $b$ in that case. At infinity it is equivalent to solving the Laplace equation with pure Neuman boundary conditions under the constraint $\int_{\Omega} u =\int_{\Omega} f$.
I would now like to have a space-variant diffusion speed $g:[0,\infty]$, with $G = \operatorname{diag}(g)$, which results in: $$d_t w(t) = -GL, \quad w(0) = b.$$ Note that this is not $d_t u=\nabla\cdot (g\nabla u)$, it is $d_t u = g\Delta u$. The solution is now $w(t) = \exp(-tGL) b$. My issue arises when $g_i \to\infty$ for some $i$. In that case I would like this to behave similar to the space-homogeneous case, i.e. at $i$ I should essentially have a Laplace equation instead. The main problem is that unlike the scalar $t$, $G$ does not commute with $V$, so I cannot do: $GV\Lambda V^T = VG\Lambda V^T$, so I cannot just plug in the minus infinity inside the exponent to get a zero.
Is there any way this can be handled in a robust manner? I am trying to implement this numerically. The only thing I came up with was multistep implicit Euler (I reparametrize $g_i$ as $g_i = \frac{c_i}{1-c_i}$, for $c_i\in [0,1]$ and I move $1-c_i$ to the left-hand side which allows me to discretise it without infinities). Since implicit Euler is uncondinitionally stable for $t\geq 0$ and $d_i\geq 0$ (since my $L$ is positive semi-definite), I can take as few steps as I want (e.g. one step). On the other hand implicit Euler with few steps poorly approximates the spectral solution even in the space-homogeneous case (and there are technically infinitely many steps to infinity, e.g. if I were to compare to explicit Euler). Is there any numerically robust way I can find the eigenvectors and eigenvalues of $GL$ as some terms $g_i$ go to infinity and others go to zero? Standard approaches cannot really handle infinite rows, so I can't just apply a standard approach to $GL$, I need something like $PDP^{-1} =\lim_{g|_J \to \infty} GL$ for a fixed set of indices $J$. My assumption here is that the infinity would show up inside of the eigenvalues of the eigendecomposition, which I do not mind, since I plan to take the exponential of that negated.