The integral2 and integral3 functions do not allow the 'ArrayValued' name-value pair, so it is necessary to ‘nest’ integral calls if you want to do that.
Try this:
%magnetic field for an off axis point due to a thick solenoid in z-x plane
l=20; %length of the solenoid
r1 = 5; %inner radius
r2 = 10; %outer radius
u0 = 1.26*(10^-3); %permeability of free space
i = 3; %current in amps
N = 800; %number of turns
C=(u0*i*N)/(2*pi*10*l);
%location of x,y,z cooridnates of P
%rho = linspace(0,22,100);
rho = 0;
h=20:1:35;
Binz = @(r,z,phi)(r.*(r-(rho.*cos(phi))))./(((rho.^2) + (r.^2) + ((h-z).^2) - (2.*r.*rho.*cos(phi))).^(3/2));
% Bz = C*integral3(Binz,r1,r2,0,20,0,2*pi); %magnetic field in z direction
Bz = integral(@(phi) integral(@(z) integral(@(r) Binz(r,z,phi), r1,r2, 'ArrayValued',1), 0,20, 'ArrayValued',1),0,2*pi, 'ArrayValued',1)
plot (Bz,h);
producing:
Bz =
29.3733 25.2358 21.3369 17.8589 14.8792 12.3925 10.3485 8.6806 7.3222 6.2144 5.3077 4.5619 3.9451 3.4317 3.0019 2.6398
and this plot:
Please check that to be certain it appears reasonable.
Just to be sure, I checked this approach with the first example in the integral3 documentation:
fun = @(x,y,z) y.*sin(x)+z.*cos(x)
q = integral(@(z) integral(@(y) integral(@(x) fun(x,y,z), 0, pi, 'ArrayValued',1), 0,1, 'ArrayValued',1), -1,1, 'ArrayValued',1)
with:
q =
2.0000
being the same as the result in the integral3 documentation.
