In the course of a discussion in the chat there emerged an interesting problem, namely a particle in an infinite well with an additional Dirac-delta function spike of scalable hight: $$ H = -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V_0\delta(x-\frac{L}{2}) $$ and with the box walls at $x=0$ and $x=L$. The solution of the problem (see linked page) ultimately pans down to the solution of the equation $$ \cot(kL/2) = -\frac{\kappa}{2k}, $$ with a characteristic wave vector $\kappa=\frac{2m}{V_0\hbar^2}$. The ground state engery would then evaluate to $$ E_0 = k_0^2\frac{\hbar^2}{2mL^2}$$ with $k_0$ the first non-trivial non-negative solution of the cotangens equation above. Interestingly the solutions of this equation are under mathematical debate. It seem to to be not clear how "badly" transcendental they are (see here and linked therein pages for more on that).
However, since the basic idea of the design of the problem was, that it should be easily accessible by perturbation theory and since that might even result in some insights into to mathematical aspects of the referred equation, the question is how to calculate the ground state energy using perturbation theory, and namely in arbitrary order.
My attempt: The obvious choice for the unperturbed system is the particle in the box, i.e. $V_0 = 0$ and to simplify things I first choose the box lengths $L$ to be $\pi$. In this way the well known state functions are
$$ \psi_k = \sqrt{\frac{2}{\pi}}sin(k x) $$
for $0 \le x \le \pi$ with the corresponding energies
$$ E_k=\frac{\hbar^2 k^2}{2m} $$ From perturbation theory we obtain a series expansion of the desired ground state energy $E_0^p$ of the perturbed system:
$$E_0^p = E_0^{(0)} + E_0^{(1)} + E_0^{(2)} + E_0^{(3)} ... $$
with 
such that we have to set $n=0$ in these expressions in order to obtain the ground state energy (note that doubly occuring indices $k_j$ in these terms imply summation over all except $j=0$, which correpsonds to $n(=0)$).
Since our perturbation (potential) is $V_0$ we have to use the following identifications:
\begin{eqnarray*}
V_{k_1 k_2} & = & <\psi_{k_1} | V_0 | \psi_{k_2}>
\\[3mm]& = &
\int_{0}^{\pi} \sqrt{\frac{2}{\pi}}\sin(k_1 x)\delta\left(x-\frac{\pi}{2}\right)\sqrt{\frac{2}{\pi}}\sin(k_2 x)dx\\[3mm]
& = &
\frac{2V_0}{\pi}\sin(k_1 \pi/2)\sin(k_2 \pi/2) \\
\end{eqnarray*}
(the latter two factors evaluate either to one of $-1,0,1$) and from the energy levels of the particle in the box
$$ E_{{k_1}{k_2}}=E_0(k_{1}^2 - k_2^2)$$
With $c=\frac{V_0}{E_0}\frac{2}{\pi}$ and evaluating the infinite sums (we note that to the quadratic eigenvalue law and the bound integrals in the enumerator convergence is no issue) in the perturbation expansion (I did that using mathematica where they all seem to evaluate nicely) one obtains:
\begin{eqnarray*} E_0^p & = & E_0\left(1+c-c^2\frac{1}{4}+c^3\left(\frac{1}{16}-\frac{\pi^2-9}{48}\right)+c^4\left(-\frac{1}{64} +\frac{\pi^2-9}{192}-\frac{\pi^2-9}{96}+\frac{\pi^2-10}{64}\right) + \mathcal{O}(c^5)\right) \\ \end{eqnarray*}
Setting here the "hight" $\frac{V_0}{E_0}$ of the spike to $\frac{1}{3}$ one obtains \begin{eqnarray*} n & & E_0^{p(n)}/E_0 \\ & & \\ 0 & & 1 \\ 1 & & 1.2122 \\ 2 & & 1.20095 \\ 3 & & 1.20137 \\ 4 & & 1.20137 \\ ex. & & \mathbf{1.20136} \\ \end{eqnarray*}
where $ex.$, the exact one is $\frac{4}{\pi}x_0^2$ from the first non-zero non-trivial numerical solution $x_0$ of $$\tan x = -\lambda x$$ with $\lambda=\frac{4E_0}{\pi V_0}$, which is equivalent to the cotangens equation from above, with $kL/2\to x$.
So that looks quite good. Now my actual question:
I would like to access arbitrarily high orders $n$ of $E_0^{(n)}$. Which would be a convenient and possibly easily implementable way to do this (e.g. using a CAS)? I have read that there are procedures to derive the terms by recursion and using derivatives (which would makes sense, since its a kind of Taylor (McLaurin) expansion). Since it seems to be a particularly simple expansion I wonder if it might even be possible to derive analytic expressions for the terms.
Here is plot of the 5 perturbation expansions, with the relative delta spike height $V_0/E_0$, that is the size of the perturbation on the $x$ axis and the relative energy on the $y$ axis ($y=1$ for the unperturbed system). The perturbation expansions are: 0-order blue, 1st-order yellow, 2nd order green, 3rd order red, 4th order purple and the "exact" in brown.
