Perturbation theory with infinite potential

It will help you to think in terms of V'(X) - V(X) .... Where is this function zero and where is it not? Where it is not, note that it takes a familiar form (up to a minus sign). It may look "bad," however your first order energy shift will involve the unperturbed wave function - for the infinite well you will find, then, that it is easy to work out the correction.....

There is a related problem that sugests how to to proceed. I learned it from Alan MacKane and I think he told me that he found it in a paper by Sidney Coleman:

Let $psi_lambda(x)$ be a normalized bound-state solution to the Schr{"o}dinger equation on the entire real line $$ [-partial_x^2 +q(x)]psi_lambda=lambda psi_lambda, quad psi_lambdain {rm L}^2({mathbb R}), $$ with $ psi_lambda(x) sim Ae^{-beta|x|}$ at large $|x|$. Then the small shift in $lambda$ that arises from confining the system in a box of length $L$, so that $psi_{lambda+deltalambda}(L/2)=psi_{lambda+deltalambda}(-L/2)=0$, is $$ delta lambda sim 4A^2 beta e^{-beta L}. $$ This is another case where the potential perturbation is infinite even though $1/L$ is a small parameter.

To obtain this result, we proceed as follows: We need to change $lambda$ so that a solution that is zero at $x=-L/2$ evolves to zero at $x=+L/2$. To do this we require a Green function for the initial value problem $$ [-partial_x^2 +q(x)-lambda]psi=f(x),qquad psi(-L/2)=0. $$ To construct the Green function we make use of the second solution
$$ chi_lambda(x) propto psi_lambda(x) int^x_0 [psi_lambda(xi)]^{-2} dxi, $$ which we scale so that $$ chi_lambda(x)sim {rm sgn}(x) A e^{beta |x|}, quad |x| gg 0. $$ The initial-value Green function is then $$ G(x,xi) = frac{1}{2A^2beta}[ psi_lambda(x) chi_lambda(xi)-chi_lambda(x) psi_lambda(xi)] theta(x-xi), $$ where $2A^2beta$ is the Wronskian $W[psi_lambda, chi_lambda]$.

Now
$$ phi(x)= frac 1{2A} left(e^{beta L/2}psi_lambda(x)+e^{-beta L/2}chi_lambda(x)right) $$ is zero at $x=-L/2$ and unity at $x=+L/2$.
We use $G(x,xi)$ with $f(x) to deltalambda, phi(x)$ to solve, to first order in $deltalambda$, the equation with $lambda to lambda+delta lambda$ and $ psi_{lambda+deltalambda}(-L/2)=0$. We obtain
$$ psi_{lambda+deltalambda}(x) = phi(x) +delta lambda int_{-L/2}^x frac{1}{2A^2beta}[ psi_lambda(x) chi_lambda(xi)-chi_lambda(x) psi_lambda(xi)] phi(xi),dxi. $$

By construction this $psi_{lambda+deltalambda}(x) $ satisfies the boundary condition at $x=-L/2$, while at $x=+L/2$ we have $$ psi_{lambda+deltalambda}(L/2) = phi(L/2)+deltalambda int_{-L/2}^{L/2}frac{1}{2A^2beta} left[ textstyle{frac 12} e^{-beta L}chi_lambda^2(xi)- textstyle{frac 12} e^{beta L} psi^2_lambda(xi)right],dxi. $$ Now $phi(L/2)=1$, $$ e^{-beta L} int_{-L/2}^{L/2} chi_lambda ^2(xi),dxi $$ is $O(1)$ and so is negligible compared to the $e^{beta L}$ term, and $$ int_{-L/2}^{L/2} psi^2_lambda(xi),dxisim 1 $$ by the normalization assumption. Thus $$ psi_{lambda+deltalambda}(L/2)= 1 - frac{deltalambda}{4A^2beta} e^{beta L}+O(deltalambda^2). $$ Requiring this to be zero gives $$ delta lambda = 4A^2 beta e^{-beta L}, $$ as claimed.

So you can do essentially the same thing for your $epsilon$ boundary shift.