Numerical integration of Airy functions

I have a problem when I try to numerically integrate the Airy function at "extreme values". Here's a MWE:

clc
clear 
%%Parameters
M=4.5
lambda_0=0.332;
tic
for pp=4:4
for m=1:1
beta1=(pp-1)/10;
alpha1=m;
omega1=7
gamma1=sqrt((M^2-1)*alpha1^2-beta1^2);
N=25;
YMAX=80;
eta01=-i*omega1/((i*alpha1*lambda_0)^(2/3))
eta_inf1=((i*alpha1*lambda_0)^(1/3))*YMAX+eta01;
%%Operations with Airy function
syms lu
aiprime(lu) = airy(1,lu);
aisecond(lu)= diff(airy(1,lu));

airy_D1=airy(1,eta01)
airy_DD1=vpa(aisecond(eta01));
airy_INT1=integral(@(n) airy(n),eta01,eta_inf1)
%%Reflection coefficient
inte1=alpha1*(alpha1^2+beta1^2)*airy_INT1;
deri1=gamma1*lambda_0*((i*alpha1*lambda_0)^(2/3))*airy_D1;
Ref_coeff1=-(inte1+deri1)/(inte1-deri1);
q1=(alpha1^2+beta1^2)*(1+Ref_coeff1)/((i*alpha1*lambda_0)^(2/3)*airy_D1);
p1=1+Ref_coeff1;
%%Coefficients
aa1=2*pi*complex(0,1)*gamma1*beta1*lambda_0*airy_D1/(alpha1*(alpha1^2+beta1^2)*airy_INT1-gamma1*lambda_0*airy_D1*(i*alpha1*lambda_0)^(2/3));
bb1=-aa1*airy(2,eta01)*airy_INT1/airy(eta01);

eta1=@(y) ((i*alpha1*lambda_0)^(1/3))*y+eta01;

Gi1=@(x) -(airy(2,x)*integral(@(n) airy(n),eta_inf1,x)-airy(x)*integral(@(n) airy(2,n),eta01,x));

V1=@(eta) ((i*alpha1*lambda_0)^(-1/3))*q1*integral2(@(x, y) airy(y),...
        eta01, eta, eta01, @(x) x);

DV1=@(eta) q1*integral(@(n) airy(n),eta01,eta);

W1=@(eta) aa1*Gi1(eta)+bb1*airy(eta);
U1=@(eta) i/alpha1*DV1(eta)-beta1/alpha1*W1(eta);

toc
end

end

ve=0:0.01:N;

U1VEC=arrayfun(U1,arrayfun(eta1,0:0.01:N));
V1VEC=arrayfun(V1,arrayfun(eta1,0:0.01:N));
W1VEC=arrayfun(W1,arrayfun(eta1,0:0.01:N));

figure(1)

subplot(1,3,1)
plot(abs(U1VEC),ve,'-k','LineWidth',2)
xlabel('$\hat{U}_1$','Interpreter','latex','Fontsize',15)
ylab=ylabel('$Y$','Interpreter','latex','Fontsize',15,'rot', 0);
set(ylab, 'Units', 'Normalized', 'Position', [-0.4, 0.45, 0]);
set(gca,'linewidth',1)
subplot(1,3,2)
plot(abs(V1VEC),ve,'-k','LineWidth',2)
xlabel('$\hat{V}_1$','Interpreter','latex','Fontsize',15)
set(gca,'linewidth',1)
subplot(1,3,3)
plot(abs(W1VEC),ve,'-k','LineWidth',2)
xlabel('$\hat{W}_1$','Interpreter','latex','Fontsize',15)
hold on
plot(abs(-i*beta1*p1/((i*alpha1*lambda_0)^(2/3))./eta1(5:0.01:N)),5:0.01:N,'--k')
set(gca,'linewidth',1)

The problem is concretely in

Gi1=@(x) -(airy(2,x)*integral(@(n) airy(n),eta_inf1,x)-airy(x)*integral(@(n) airy(2,n),eta01,x));

In this case, eta01 and eta_inf1 take very extreme values, making the U1 and W1 solutions get wrong values.

Do you have any idea on how to solve this problem? Without increasing the computational time too much, so avoiding increasing the precision artificially with vpa.

[John BG]

Hi Carlos

this is John BG [e-mail]

in the definition of anonymous function Gi1, why do you define some airy functions as expected, with 2 inputs

airy(2,x)

yet you only use one argument in some other appearances of same airy function while defining Gi1?

the default type is 0, so if you miss the 1st argument MATLAB will automatically use

airy(0,x)

wherever the function type index is omitted.

[carlos g]

Hi John BG,

Yes, I did it on purpose. Sometimes I am using the Ai function and, in the first case you mentioned, it's the Bi function. Is this your question?

[carlos g]

Any idea? Thank you