code 35154;
  procedure NONEXP SPHER BESS I(X, N, I); value X, N;
  real X; integer N; array I;
  if X= 0 then 
  begin I[0]:=1;
    for N:= N step -1 until 1 do I[N]:= 0
  end else 
  begin real X2, R, S; integer M;
    X2:= X+X;
    I[0]:= X2:= if X = 0 then 1 else if X2 < 0.7 then 
        SINH(X) / (X * EXP(X)) else (1-EXP(-X2))/X2;
    if N= 0 then "GO TO" EXIT;
    R:= 0; M:= START(X,N,1);
    for M:= M step -1 until 1 do 
    begin R:= 1/((M+M+1)/X+R);
      if M  <=  N then I[M]:= R
    end;
    for M:= 1 step 1 until N do 
    I[M]:= X2:= X2 * I[M];       EXIT:
  end;
        eop