Try this instead:
ALP = @(j,l,x) ((-1)^l)*(1-x*x)^(l/2)*diff(legendreP(j,x),l);
syms x y
alpha = 45*pi/180;
b = cos(alpha);
l = 0;
N = 10;
phi = 0;
X = (sqrt(1-x*x)/sin(alpha))*cos(phi);
Y = (sqrt(1-x*x)/sin(alpha))*sin(phi);
Z = (x-cos(alpha))/sin(alpha);
r = sqrt(X*X+Y*Y+Z*Z);
y1 = x;
y2 = (x/2)*log((1+x)/(1-x)) - 1;
for i=1:N
for j=1:N
fun = y1.*int(subs(y2.*r^j,x,y).*ALP(j,l,y),y,x,1);
G0(i,j) = vpaintegral(fun,x,b,1)
end
end
The integral function will not work with symbolic variables or symbolic functions. The vpaintegral function (R2016b and later releases) will. Use the double function with ‘G0’ to convert it to a double array (with some loss of precision) if you want to use it in non-symbolic calculations.
