% compute the first i eigenfunction of stable spline kernel 
% in the interval [0,T] sampled with interval dT
% beta is the decay rate (exponent of change of coordinates 
% tau = exp(-beta*t)
% the j-th column of Phi is the j-th eigenfunction
% lambda are the eigenvalues

function [Phi lambda] = stable_spline(i,beta,T,dT); 

t = [0:dT:T];

for k=1:i,
    
%alpha = solve('1/cosh(x-(num2str(k)-pi/2))+cos(x-(num2str(k)-pi/2))=0','x');
alpha = fzero(@(x) fun_eig(x),(k-1/2)*pi);
lambda(k)=(1/alpha)^4;
C31=(2/(1+exp(-2*alpha))-1)/sin(alpha);
C21=1-C31*exp(-alpha);
C11=-1-C31*exp(-alpha);

dtau=0.00001;
tau=[0:dtau:1];
f_int = C11*cos(alpha*tau)+C21*sin(alpha*tau) + C31*exp(-alpha*(1-tau))+exp(-alpha*tau);
f_int=f_int.^2;
C4 = 1/sqrt(sum(f_int)*dtau);
C3=C4*(2/(1+exp(-2*alpha))-1)/sin(alpha);
C2=C4-C3*exp(-alpha);
C1=-C4-C3*exp(-alpha);

Phi(:,k)=C1*cos(alpha*exp(-beta*t))+C2*sin(alpha*exp(-beta*t))+C3*exp(-alpha*(1-exp(-beta*t)))+C4*exp(-alpha*exp(-beta*t));
end;
%Phi = Phi(end:-1:1,:);

