function [keps,klam] = heunpt(n0, m, N)
% Returns the perturbative expansion of any eigenvalue of the 
% particle on the singular surface defined in Fateev-Onofri, july2003
% (FermilabPub-03/165-T, LPM/03-013, UPRF-2003-12)
% An iterative method is used to calculate the perturbative expansion of
% the "accessory parameter"; this is transformed into a series for kappa
% and also a new series in the parameter lambda=(1+exp(-4u))^(-1).
% Usage: [keps,klam] = heunpt(n, m, N)
%     N=pert-order (unfeasible to go beyond 500)
%     keps=coeff for \kappa(u)\approx c_n eps^n
%     klam=coeff for \kappa(u)\approx c_n lambda^n
% Note: Asymmetric formulation  (unnormalized Jacobi polynomials)
% © E.Onofri and V.Fateev, Fermilab, June 2003, Parma July 2003

warning off MATLAB:divideByZero

Nmax = N+n0+1;

delta=zeros(N+1,1);
delta(1) = n0*(n0+m+1);
n = (1:Nmax)';
ns=(0:Nmax)';

ad  = spdiags(2*n.^2.*(m+n).^2./((m+2*n-1).*(m+2*n)),-1,Nmax,Nmax);
a   = spdiags(2*(n-1).^2.*(m+n-1).^2./((m+2*n-2).*(m+2*n-1)), 1,Nmax,Nmax);    

R      =  (ns-n0).*(ns+n0+m+1);
R(n0+1)=  1.0;
R      =  1.0./R;
R(n0+1)=  0.0;

R  = spdiags(R, 0, Nmax, Nmax);   % risolvente nel
                                  % sottospazio ortogonale 

if m>0,
  Vd = m*(m*(1+m)+(1+m)*(2+m)*ns+(2+m)*ns.^2)./((m+2*ns).*(m+2*ns+2)); 
else
  Vd = zeros(size(ns));
end

V = a+ad+spdiags(Vd,0,Nmax,Nmax);	% perturbation off--diagonal + diagonal

v{1} = sparse(Nmax,1); 
v{1}(n0+1) = 1; % questo è |n>/<n|n>
eta{1} = v{1};
w0 = v{1}'*V; 

delta(2) =  w0*eta{1};

for(j=2:N),				% start iteration
  
  Wtmp = -V*eta{j-1};
  
  for(k=1:j-1),	
		
    Wtmp = Wtmp+delta(k+1)*eta{j-k};
  
  end

  eta{j}     = R*Wtmp;
  delta(j+1) = w0*eta{j};

  if(isnan(delta(j+1))),
    break;
  end
  
end

%% keps = (48qx-6(2m+1)(3x-1))/(1+x) --> next three lines 
keps = 24*delta;
keps(1)=keps(1)+6*(1+2*m);
keps(2)=keps(2)+6*(1-2*m);

keps=invert(keps);

klam = [keps(1);convert(keps(2:end))];	% skip first coeffcient

%--------------------------------
% Invert (1+x) K = Y(x)
function S=invert(Y)

N=length(Y);
Xp1 = speye(N)+spdiags(ones(N,1),-1,N,N);
S=Xp1\Y;

% To manipulate power series, use sparse matrix representation:
% X=spdiags(ones(N,1),-1,N,N);
% e.g. (1+X)K=D -> K=(1+X)\D (left matrix division=mldivide)

%---------------------------------

function C = convert(K)
% transforms the series in epsilon into the series in lambda:
% epsilon=lambda/(2-lambda). The relation is governed by Pascal 
% scheme of binomial coefficients.

N=length(K);
n=(1:N)';
Norm=spdiags(2.^(-n),0,N,N);
S = Norm*sparse(abs(pascal(N,1)));
C = S*K;
