X=-2:0.1:6;
y=0.4;
t=0.75;
v=1.5;
alpha=0.1;
for i=1:numel(X)
    x = X(i);
    Fun = @(r)fun(r,x,y,t,v,alpha);
    C(i) = integral(Fun,0,Inf);
end
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   5.9e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   6.5e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   6.7e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   7.1e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   7.4e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   8.1e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   8.3e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   9.4e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   9.7e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   1.0e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   1.1e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   1.2e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   1.2e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   1.3e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   1.4e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   1.6e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   1.6e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   1.7e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   1.9e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   2.0e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   2.2e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   2.3e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   2.4e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   2.6e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   2.8e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   3.0e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   3.2e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
plot(X,C)

function value = fun(r,x,y,t,v,alpha)
  gamma=1/(1-alpha);
  A=exp(v/2*(x+y)-alpha*gamma*t);
  B=(r.^3)./(r.^2+gamma+v^2/2).^2;
  D=exp((alpha*gamma.^2*t)./(r.^2+gamma+v^2/2));
  K=besselj(0,r.*sqrt(x.^2+y.^2));
  L=besselj(2,r.*sqrt(x.^2+y.^2));
  M=B.*D.*(K+L);
  value = alpha*gamma^2*y*A*M;
end




