You have
Bz= @(x,y) (((u./(4.*pi)).*((((((-1).^1).*d1(ai))./(r1(ai).*(r1(ai)+((-1).^(1+1)).*C1(ai))))-(C1(ai)./(r1(ai).*(r1(ai)+d1(ai)))))+(((((-1).^2).*d2(ai))./(r2(ai).*(r2(ai)+(((-1).^(2+1))).*C2(ai))))-(C2(ai)./(r2(ai).*(r2(ai)+d2(ai)))))+(((((-1).^3).*d3(ai))./(r3(ai).*(r3(ai)+(((-1).^(3+1))).*C3(ai))))-(C3(ai)./(r3(ai).*(r3(ai)+d3(ai)))))+(((((-1).^4).*d4(ai))./(r4(ai).*(r4(ai)+((-1).^(4+1)).*C4(ai))))-(C4(ai)./(r4(ai).*(r4(ai)+d4(ai))))))))
The calculations on the right hand side do not use x or y, so whatever it returns is not going to be the same size as x and y
The function you pass into integral2() is not passed scalar x and y: it is passed arrays of x and y that are the same size, and it is responsible for returning a result that is the same size as x.
