function[F]=mbssel(X1)

% The function MBSSEL generates a zero order modified Bessel function of the first kind:
%     Io(X) = SUM( (X/2)**(2*j)/(j!)**2 )  j=0,1,2,3,...
%
% Calling sequence-
% [F]=mbssel(X1)
% 
% Input-
%	x1	- a square matrix x1(n,n)
% Output-
%	f	- a value of a Bessel function
%
% See also-
% in ORMSBY it is used as [f0]=mbssel(5).

%----- Generate a function
F=1.+(X1^2/4) ;

for  M=2:10       ;
   J=1            ;
   M1=M-1         ;      
   for I=1:M1     ;  % get j! 
       J=J*(I+1)  ;
   end   
   RJ=J           ;
   F=F+(X1^(2*M))/(2^(2*M)*RJ^2);  
end