function J2=jacobian_map_phi_even(beven)
global mu sigma
%aodd=2*bodd;
m=length(beven);
c2fft=quadratic_sumFFT_real([beven;zeros(m,1)]);

% tic
% J2=zeros(m,m);
% 
% for n=1:m
%     for j=1:n
%             J2(n,j)=(pi^2*(2*n)^2+sigma*mu)*(n==j)+real(6*sigma*(c2fft(n-j+1)-c2fft(n+j+1)));
%             J2(j,n)=J2(n,j); 
%     end
% end
% toc
%tic
J2bis=zeros(m,m);
for n=1:m
    J=1:n;
    J2bis(n,J)=real(6*sigma*(c2fft(n-J+1)-c2fft(n+J+1)));
end
Diag=diag(pi^2*(2*[1:m]).^2+sigma*mu);
J2bis=J2bis+Diag;
J2=J2bis+J2bis.'-diag(diag(J2bis));

    


end

function [s]=quadratic_sumFFT_real(a1)

m=length(a1);

b1=-flipud(a1(1:m));
ta1=[zeros(m,1);b1;0;a1;zeros(m,1)];
tu1=ifft(ifftshift(ta1));

Ft=fftshift(fft(tu1.*tu1));

s=(4*m+1)*Ft(2*m+1:3*m+1);
end

function [s]=quadratic_sumFFT_real_int(a1)

i2=intval('2');
i4=intval('4');

a=[intval(0);a1];
M=size(a,1);
% We make sure that the inputs are powers of 2 %
M1=(2^ceil(log(4*M-1)/log(2))-4*M)/2;   % Hence 4*M+2*M1 is a power of 2 %
ta1=[intval(zeros(M+M1+1,1));-a(M:-1:2,1);a;intval(zeros(M+M1,1))];
tu1=verifyfft(ifftshift(ta1),-1);
F=fftshift(verifyfft(tu1.*tu1,1));
s=(i4*intval(M)+i2*intval(M1))*F(M1+2*M+1:M1+3*M,1);
end

% function [s]=quadratic_sum_loops(a1)
% 
% m=length(a1);
% a1=[a1;0];
% s=zeros(m,1);
% 
% for k=0:m-1  
%     s(k+1,1)=0;
%     for l=-(m-1):m-1
% 
%         if ( abs(k-l) <= m ), 
%             s(k+1,1)=coeff(a1,l)*coeff(a1,k-l)+s(k+1,1);
%         end
%         
%     end
% end
% 
% end
% function [coef]=coeff(c,n1)
% if n1<0
%     coef=-(c(abs(n1)+1));
% else
%     coef=c(abs(n1)+1);
% end
% end


