算个东西，需要用到bessel函数的根。找到了一个maple的算法，然后改成matlab的。如下：
maple代码的页面：http://www.math.rwth-aachen.de/mapleAnswers/html/155.html
代码：
> ZerosBesselJ := proc (maxv, maxs)
> local j, incr, v, h, s;
>
> j := array(0..maxv, 1..maxs);
> incr := 4.0;
> for v from 0 to maxv do
> h := evalf(v + 1.9*v^(1/3) + 1);
> if v = 0 then
> j[v,1] := fsolve(BesselJ(v,x), x, 2.0 .. 3.0)
> else
> j[v,1] := fsolve(BesselJ(v,x), x, 2.0 .. h)
> fi;
> for s from 2 to maxs do
> j[v,s] := fsolve(BesselJ(v,x), x, j[v,s-1] .. j[v,s-1]+incr)
> od
> od;
> RETURN( eval(j) )
> end:
改成matlab之后的代码如下：
clear all;
maxv = 10;
maxs = 10;
j= zeros(maxv, maxs);
incr = 4.0;
for v=0:maxv-1
h = v+1.9*v^(1/3)+1;
if (v==0)
j(v+1,1) = fzero(@(x)besselj(v,x),2);
else
j(v+1,1) = fzero(@(x)besselj(v,x),h);
end
for s=2:maxs
j(v+1,s) = fzero(@(x)besselj(v,x),j(v+1,s-1)+incr);
end
end
j
运算结果：
j =
2.4048 5.5201 8.6537 11.7915 14.9309 18.0711 21.2116 24.3525 27.4935 30.6346
3.8317 7.0156 10.1735 13.3237 16.4706 19.6159 22.7601 25.9037 29.0468 32.1897
5.1356 8.4172 11.6198 14.7960 17.9598 21.1170 24.2701 27.4206 30.5692 33.7165
6.3802 9.7610 13.0152 16.2235 19.4094 22.5827 25.7482 28.9084 32.0649 35.2187
7.5883 11.0647 14.3725 17.6160 20.8269 24.0190 27.1991 30.3710 33.5371 36.6990
8.7715 12.3386 15.7002 18.9801 22.2178 25.4303 28.6266 31.8117 34.9888 38.1599
9.9361 13.5893 17.0038 20.3208 23.5861 26.8202 30.0337 33.2330 36.4220 39.6032
11.0864 14.8213 18.2876 21.6415 24.9349 28.1912 31.4228 34.6371 37.8387 41.0308
12.2251 16.0378 19.5545 22.9452 26.2668 29.5457 32.7958 36.0256 39.2404 42.4439
13.3543 17.2412 20.8070 24.2339 27.5837 30.8854 34.1544 37.4001 40.6286 43.8438