I don't know what bessel0 is, so I just substituted atan() below, but this the way you would do that with numerical integration. Notice that final() is now a function of 3 variables. It is easiest to first define a function that will only accept a scalar z value and then use arrayfun to make it work with arrays of z values.
m1=0.00098;
m2=0.91426;
lambidazero=659.57;
sigma=10.981;
deltan=-0.59;
r=1;
% x from 500 until 750 and ro from 0 until r
final = @(x,ro,z)(m1 + m2.*exp(-((x-lambidazero).^2)./((sigma).^2))) ...
.* (exp((-1i.*deltan.*2.*pi.*(((r.^2)-(ro.^2)).^(1./2)))./(x))) ...
.* atan(z.*ro).*ro;
fscalar = @(z)integral2(@(x,y)final(x,y,z),500,750,0,r);
% The function fscalar(z) can only accept a scalar value for z.
% Use arrayfun to make a function that accepts vectors.
f = @(z)arrayfun(fscalar,z);
% Now plot f just as you would any function in MATLAB
x = 0:0.1:1;
y = f(x);
plot(x,y)
