function m550x

% Compute integral on a representative grid of k-values:
kk=[]; II=[];
for nk=1:99;
    k=(nk-1)/100;
    II(nk)=integral(@(x) 1./abs(1-k*x.^3).^(3/2),0,1);
    kk=[kk k];
end

plot(kk,II,'k*')
xlabel('k')
ylabel('I(k)')
grid
hold on

%  Now for the approximations:
a0 = 1; a1 = 3/8; a2 = 15/56; a3 = 7/32;
plot(kk,a0+a1*kk+a2*kk.^2+a3*kk.^3,'k--')    % Small k
plot(kk,2/3./(1-kk).^(1/2),'r--')            % k -> 1

res = [kk(91) II(91)];

res = [res a0+a1*kk(91)+a2*kk(91).^2+a3*kk(91).^3 2/3./(1-kk(91)).^(1/2)];

%  Multiplicative extraction:
%    I(k) = a_0 + a_1 k + a_2 k^2 + a_3 k^3
%         = (1-k)^{-1/2} (g_0 + g_1 k + g_2 k^2 + g_3 k^3)
%  which gives
g0 = a0; g1 = a1-a0/2; g2 = a2-a1/2-a0/8; g3 = a3-a1/8-a2/2-a0/16;
plot(kk,(g0+g1*kk+g2*kk.^2+g3*kk.^3)./(1-kk).^(1/2),'b-')
% limit for k -> 1 is AA*(1-k)^{-1/2} with
AA=g0+g1+g2+g3;
% versus
A=2/3;
[AA A]
% (so the approach to the singular limit is not quite correct)

res = [res (g0+g1*kk(91)+g2*kk(91)^2+g3*kk(91)^3)/(1-kk(91))^(1/2)];

%  Additive extraction:
%    I(k) = A (1-k)^{-1/2} + g_0 + g_1 k + g_2 k^2 + g_3 k^3
%  which gives
g0 = a0-2/3; g1 = a1-1/3; g2 = a2-1/4; g3 = a3-5/24;
plot(kk,2/3./(1-kk).^(1/2)+g0+g1*kk+g2*kk.^2+g3*kk.^3,'r-')

res = [res 2/3/(1-kk(91))^(1/2)+g0+g1*kk(91)+g2*kk(91)^2+g3*kk(91)^3]

hold off
axis([0 1 .95 4])

  h=legend('Numerical','$k\to0$','$k\to1$','Multiplicative','Additive',...
	   'Location','northwest');
set(h,'interpreter','latex','fontsize',12)

end
