input_case nil$

m:=6$

% Quantum nonlinear oscillator
% ----------------------------
% perturbation theory to the order m
array d(m); operator x,c,de; % x(k,j) = <n+j|x^k|n>
x(1,1):=sqrt((n+1)/2)$ x(1,-1):=sqrt(n/2)$ x(1,0):=0$
for k:=2:m+2 do for j:=-k:k do x(k,j):=
    (if j<k-1 then sub(n=n+j,x(1,1))*x(k-1,j+1) else 0)
    +(if j>1-k then sub(n=n+j,x(1,-1))*x(k-1,j-1) else 0);
for k:=3:m+2 do for j:=-k:k do x(k,j):=c(k-2)*x(k,j);

procedure v(l,j,a); % perturbation theory
<<  if l<m then
    for k:=2 step 2 until fix((3*l-abs(j))/3) do
        v(l-k,j,a*de(k)/j);
    for k:=1:min(l,1+(3*l-abs(j))/2) do
    for i:=max(-2-k,-3*(l-k)-j) step 2 until min(2+k,3*(l-k)-j) do
    if j+i=0
    then d(m-l+k):=d(m-l+k)+a*sub(n=n+j,x(k+2,i))
    else v(l-k,j+i,-a*sub(n=n+j,x(k+2,i))/(j+i))
>>$

v(m,0,1)$

factor c;
for k:=2 step 2 until m do
<< de(k):=d(k); write factorize(de(k)); >>;

showtime;

off nat,pri,echo; out "quant.res";
for k:=2 step 2 until m do write de(k):=de(k);
write "end"; shut "quant.res";

end;