function J = BesselFunction(order,m_max)

% 

J = zeros(1,2 * m_max + 1 + order);
Jtemp = zeros(1,2 * m_max + 1);
for m = 0:m_max
    Jtemp((2 * m_max + 1) - (2 * m)) = ((-1) ^ m) / (factorial(m) * gamma(m + order + 1) * (4 ^ m));
end
J(1:(2 * m_max + 1)) = Jtemp;
J = ((1 / 2) ^ order)* J;

end