Bessel Matlab

.m

n = input('Digite n (El número funciones de Bessel a encontrar: ');
Bmax = input('Digite Bmax (El valor máximo de B): ');

disp('Funciones de Bessel (Usando el método de Simpson)');

fprintf('B   \t\t');
arrayfun(@(x)fprintf('J%04d\t\t', x), 0:n, 'UniformOutput', false);
fprintf('\n');

syms k v t
fJ(k, v, t) = cos(k * t - v * sin(t));

for i = 0:n
    for b = 0:Bmax * 10
        fa = double(fJ(i, b / 10, 0));
        fb = double(fJ(i, b / 10, pi));
        fhalfab = double(fJ(i, b / 10, pi / 2));
        coef = 1 / 6;
        J(i + 1, b + 1) = coef * (fa + 4 * fhalfab + fb);
    end
end

for b = 0:Bmax * 10
    fprintf('%.2f\t\t', b / 10);
    for i = 0:n
        fprintf('%.4f\t\t', J(i + 1, b + 1));
    end
    fprintf('\n');
end