function diskmodes(m,n)
% DISKMODES - creates plots of the solutions of helmholtz equation on
% a disk up to the first m=10, n=10 modes
%   
%
% check for limits on m and n
%
if (m>10 | n> 10) 
  error('Sorry: use m,n<10');
end
%
%  define the coordinate systems
%
dtheta=pi/180.; % 1 degree increments
theta=-pi:dtheta:pi;
dr=.01;
r=0:dr:1;
[R,Theta]=meshgrid(r,theta); % create matrix of (r,\theta values) see help meshgrid
%%%%%%%%%%%%%%%%%%%%5
%  now find the nth zero of Jm(z) at the moment let's do this the hard
%  way starting with the first zero and iterating forward
%%%%%%%%%%%%%%%%%%%%%%
Jm=inline(sprintf('besselj(%d,z)',m)); % make inline function for Jm(z)
if m>0 % asymptotic approximation for first zero of Jm for large m
  zm1=m+1.85575*m^(1/3)+1.03315*m^(-1/3)-.00397/m;  
else
  zm1=2.4
end
zguess=zm1 ; % estimate of first zero of J0
for i=1:n
  zmn(i)=fzero(Jm,zguess);
  zguess=zmn(i)+3;
end
str=sprintf('z_%d,%d=%f',m,n,zmn(n));
disp(str);
modemn=besselj(m,R*zmn(n)).*cos(m*Theta);
X=R.*cos(Theta); 
Y=R.*sin(Theta);
surfc(X,Y,modemn); shading interp; axis image; set(gca,'fontweight','bold','fontsize',[14]);
view(0,90);
title(sprintf('Disk Modes m=%d, n=%d,z_{mn}=%f',m,n,zmn(n)));
return