I have considered the Sturm–Liouville problem in a form of one-dimensional Schrodinger equation with Dirichlet boundary conditions: $$ \left\{\begin{array}{l} \displaystyle-\frac{\hbar^2}{2m}\varphi''(x)+V(x)\varphi(x)=E\varphi(x),\quad a\leqslant x\leqslant b,\\ \varphi(a)=\varphi(b)=0.\end{array}\right. $$ If we obtained some specific eigenvalue $E_n$ and eigenfunction $\varphi_n(x)$, then the eigenvalue $E_n$ should somehow depend on the boundary point locations $a$ and $b$. I have found that $$ \frac{\partial E_n}{\partial a}=\frac{\hbar^2}{2m}|\varphi'_n(a)|^2,\qquad \frac{\partial E_n}{\partial b}=-\frac{\hbar^2}{2m}|\varphi'_n(b)|^2. $$
I suppose this is a well-known result in the theory of differential equations or in the functional analysis, so could anyone give any additional information about this problem (any books or papers, who first obtained this result, what is its higher-dimensional generalizations etc.).
I am not so sure about the references you requested. Yet the generalization of your result to higher-dimensional cases is somewhat straightforward.
Consider the following stationary Schrodinger problem \begin{align} -\kappa\Delta\varphi+U\varphi&=E\varphi&&\text{in }\Omega,\\ \varphi&=0&&\text{on }\partial\Omega, \end{align} where $$ \kappa=\frac{\hbar^2}{2m}, $$ and $\Omega\subseteq\mathbb{R}^n$ is an $n$-dimensional simply connected domain with piecewise smooth boundary $\partial\Omega$. In addition, $\varphi$, if with more than one component, is by default as a column vector.
Note that the domain $\Omega$ can be characterized by a level set function $\phi=\phi(\mathbf{x})$, i.e., \begin{align} \phi&>0&&\text{in }\Omega,\\ \phi&=0&&\text{on }\partial\Omega. \end{align} Thus when considering the continuous deformation of $\Omega$, it suffices to focus on the change of $\phi$.
Using this notation, the stationary Schrodinger problem witnesses an equivalent statement:
In this way, the energy eigenvalue $E$ serves as a Lagrangian multiplier to conserve the probability. Obviously, this $E$ depends on $\phi$. Plus, since $E\in\mathbb{R}$, it is a functional of $\phi$. It can be determined as follows.
Suppose $\varphi$ is a solution to the above stationary point problem with some fixed $\phi$. Then $\varphi$ satisfies, for some $E\in\mathbb{R}$, \begin{align} -\kappa\Delta\varphi+U\varphi&=E\varphi&&\text{in }\left\{\phi>0\right\},\\ \varphi&=0&&\text{on }\left\{\phi=0\right\},\\ \int_{\left\{\phi>0\right\}}\left|\varphi\right|^2{\rm d}V&=1. \end{align} Left-multiply $\varphi^{\dagger}$ on both-sides of the governing equation, and the integration by parts yields \begin{align} E&=E\int_{\left\{\phi>0\right\}}\left|\varphi\right|^2{\rm d}V=\int_{\left\{\phi>0\right\}}\varphi^{\dagger}\left(E\varphi\right){\rm d}V\\ &=\int_{\left\{\phi>0\right\}}\varphi^{\dagger}\left(-\kappa\Delta\varphi+U\varphi\right){\rm d}V=\int_{\left\{\phi>0\right\}}\left(\kappa\left\|\nabla\varphi\right\|^2+U\left|\varphi\right|^2\right){\rm d}V. \end{align} Therefore, $$ E[\phi]=J[\varphi[\phi],\phi], $$ where $\varphi[\phi]$ means that $\varphi$, as a solution to the stationary point problem, depends on the choice of $\phi$.
To examine how the variation of $\phi$ has impact on the change of $E$, it suffices to figure out $$ \delta E[\phi]=\delta J[\varphi[\phi],\phi]. $$ Due to the difficulty in calculating this variation, consider a family of feasible $\phi$, denoted as $$ \phi=\phi(\mathbf{x},\tau) $$ with $\tau\in\mathbb{R}$, i.e., $\phi=\phi(\cdot,\tau)$ is a level set function for each fixed $\tau$, and it switches between different level set functions continuously as $\tau$ changes. With this setting, it suffices to determine $$ \frac{\rm d}{{\rm d}\tau}E[\phi(\cdot,\tau)]=\frac{\rm d}{{\rm d}\tau}J[\varphi(\cdot,\tau),\phi(\cdot,\tau)]. $$ Here $$ \varphi(\mathbf{x},\tau)=\varphi[\phi(\cdot,\tau)](\mathbf{x}) $$ stands for the solution to the stationary point problem as the level set function reads $\phi=\phi(\cdot,\tau)$.
Thanks to the Leibniz integral rule (since $\left\{\phi>0\right\}$ is also $\tau$-dependent), $$ \frac{\rm d}{{\rm d}\tau}E=\int_{\left\{\phi>0\right\}}\frac{\partial}{\partial\tau}\left(\kappa\left\|\nabla\varphi\right\|^2+U\left|\varphi\right|^2\right){\rm d}V+\int_{\left\{\phi=0\right\}}\left(\kappa\left\|\nabla\varphi\right\|^2+U\left|\varphi\right|^2\right)\mathbf{v}\cdot{\rm d}\mathbf{S}, $$ where $\mathbf{v}$ is the normal "velocity" of the moving boundary $\left\{\phi(\tau)=0\right\}$ as $\tau$ changes.
The derivation henceforth would adopt these notations: \begin{align} \dot{\phi}&=\frac{\partial}{\partial\tau}\phi(\cdot,\tau),\\ \dot{\varphi}&=\frac{\partial}{\partial\tau}\varphi(\cdot,\tau). \end{align} Since $\varphi=\varphi(\cdot,\tau)$ is the solution to the stationary point problem, it satisfies \begin{align} -\kappa\Delta\varphi(\cdot,\tau)+U\varphi(\cdot,\tau)&=E[\phi(\cdot,\tau)]\varphi(\cdot,\tau)&&\text{in }\left\{\phi(\cdot,\tau)>0\right\},\\ \varphi(\cdot,\tau)&=0&&\text{on }\left\{\phi(\cdot,\tau)=0\right\}. \end{align} Specifically, the second condition implies that \begin{align} \dot{\varphi}&=-\mathbf{v}\cdot\nabla\varphi&&\text{on }\left\{\phi=0\right\}, \end{align} because any front-tracking particle $\mathbf{x}=\mathbf{x}(\tau)\in\left\{\phi(\cdot,\tau)=0\right\}$ with $\dot{\mathbf{x}}(\tau)=\mathbf{v}(\mathbf{x}(\tau),\tau)$ complies with $\varphi(\mathbf{x}(\tau),\tau)=0$, for which $$ \frac{\rm d}{{\rm d}\tau}\varphi(\mathbf{x}(\tau),\tau)=0\iff\dot{\varphi}(\mathbf{x}(\tau),\tau)+\dot{\mathbf{x}}(\tau)\cdot\nabla\varphi(\mathbf{x}(\tau),\tau)=0. $$ Thanks to these facts, \begin{align} &\int_{\left\{\phi>0\right\}}\frac{\partial}{\partial\tau}\left(\kappa\left\|\nabla\varphi\right\|^2+U\left|\varphi\right|^2\right){\rm d}V\\ &=\int_{\left\{\phi>0\right\}}\left(\kappa\nabla\dot{\varphi}^{\dagger}\nabla\phi+U\dot{\varphi}^{\dagger}\varphi\right){\rm d}V+\text{c.c.}\\ &=\left[\kappa\int_{\left\{\phi=0\right\}}\dot{\varphi}^{\dagger}\nabla\varphi\cdot{\rm d}\mathbf{S}+\int_{\left\{\phi>0\right\}}\dot{\varphi}^{\dagger}\left(-\kappa\Delta\varphi+U\varphi\right){\rm d}V\right]+\text{c.c.}\\ &=\left[-\kappa\int_{\left\{\phi=0\right\}}\left(\mathbf{v}\cdot\nabla\varphi^{\dagger}\right)\nabla\varphi\cdot{\rm d}\mathbf{S}+\int_{\left\{\phi>0\right\}}\dot{\varphi}^{\dagger}\left(E\varphi\right){\rm d}V\right]+\text{c.c.}\\ &=\left[-\kappa\int_{\left\{\phi=0\right\}}\left\|\nabla\varphi\right\|^2\mathbf{v}\cdot{\rm d}\mathbf{S}+\int_{\left\{\phi>0\right\}}\dot{\varphi}^{\dagger}\left(E\varphi\right){\rm d}V\right]+\text{c.c.}\\ &=-2\kappa\int_{\left\{\phi=0\right\}}\left\|\nabla\varphi\right\|^2\mathbf{v}\cdot{\rm d}\mathbf{S}+E\int_{\left\{\phi>0\right\}}\frac{\partial}{\partial\tau}\left|\varphi\right|^2{\rm d}V, \end{align} where $\text{c.c.}$ stands for the complex conjugate of its previous term. Further, thanks to the conservation of probability, $$ 0=\frac{\rm d}{{\rm d}\tau}\int_{\left\{\phi>0\right\}}\left|\varphi\right|^2{\rm d}V=\int_{\left\{\phi>0\right\}}\frac{\partial}{\partial\tau}\left|\varphi\right|^2{\rm d}V+\int_{\left\{\phi=0\right\}}\left|\varphi\right|^2\mathbf{v}\cdot{\rm d}\mathbf{S}. $$ Consequently, $$ \int_{\left\{\phi>0\right\}}\frac{\rm d}{{\rm d}\tau}\left(\kappa\left\|\nabla\varphi\right\|^2+U\left|\varphi\right|^2\right){\rm d}V=-2\kappa\int_{\left\{\phi=0\right\}}\left\|\nabla\varphi\right\|^2\mathbf{v}\cdot{\rm d}\mathbf{S}-E\int_{\left\{\phi=0\right\}}\left|\varphi\right|^2\mathbf{v}\cdot{\rm d}\mathbf{S}. $$ This result immediately gives $$ \frac{\rm d}{{\rm d}\tau}E=\int_{\left\{\phi=0\right\}}\left[-\kappa\left\|\nabla\varphi\right\|^2+\left(U-E\right)\left|\varphi\right|^2\right]\mathbf{v}\cdot{\rm d}\mathbf{S}. $$ Taking the boundary condition \begin{align} \varphi&=0&&\text{on }\left\{\phi=0\right\} \end{align} into consideration, the above expression eventually reduces to
This is the higher-dimensional generalization of your one-dimensional result, where $\varphi$ could have more than one component.
As an example, consider the one-dimensional case, where $\left\{\phi(\cdot,\tau)>0\right\}=\left(a+\tau,b\right)\subseteq\mathbb{R}$. In this case, $$ \left\{\phi(\cdot,\tau)=0\right\}=\left\{a+\tau,b\right\}, $$ and \begin{align} \mathbf{v}&=\frac{\rm d}{{\rm d}\tau}\left(a+\tau\right)=1&&\text{at }x=a+\tau,\\ \mathbf{v}&=\frac{\rm d}{{\rm d}\tau}b=0&&\text{at }x=b. \end{align} Further, the outward unit normal of $\left(a+\tau,b\right)$ reads \begin{align} \mathbf{n}&=-1&&\text{at }x=a+\tau,\\ \mathbf{n}&=1&&\text{at }x=b. \end{align} Thus when $\tau=0$, $$ \frac{\rm d}{{\rm d}\tau}E=-\kappa\int_{\left\{\phi=0\right\}}\left\|\nabla\varphi\right\|^2\mathbf{v}\cdot{\rm d}\mathbf{S}=-\kappa\int_{\left\{\phi=0\right\}}\left|\varphi'(x)\right|^2\mathbf{v}\cdot\mathbf{n}{\rm d}x=\kappa\left|\varphi'(a+\tau)\right|^2=\kappa\left|\varphi'(a)\right|^2, $$ which tells the energy variation upon the change of the left end-point of the domain $\left(a,b\right)$. This is in agreement with your result. Likewise, set $\left\{\phi(\cdot,\tau)=0\right\}=\left(a,b+\tau\right)$, and you could figure out the energy variation upon the change of the right end-point of $\left(a,b\right)$, which also agrees with your result.