The mathematical reason is that when m = 1:2, the computation of j.*cos(m*pi)-j_minus is effectively equal to zero—absolute zero!
%Neumann Bessel Function
x=(0:0.01:15).';
i=0:50;
hold on
set(gca,'YLim',[-1,1]); %range
set(gca,'XLim',[0,15]); %domain
for m=0:2
j=sum((((-1).^i)./(factorial(i).*factorial(i+m))).*(x./2).^(2*i+m),2);
if m == 0
y(:,m+1)= bessely(m, x);
else
i_minus = m:50;
j_minus = sum((((-1).^i_minus)./(factorial(i_minus).*factorial(i_minus-m))).*(x./2).^(2*i_minus-m),2)
Term1 = j.*cos(m*pi)
Term2 = sin(m*pi)
Term3 = Term1 - j_minus
y(:,m+1)= (Term1 - j_minus)./Term2;
end
plot(x, y(:,m+1))
end
legend('Y0','Y1','Y2')
