算个东西，需要用到bessel函数的根。找到了一个maple的算法，然后改成matlab的。如下：

代码：

> 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