function isosurf(x,y,z,u,v,c)
%isosurf(x,y,z,U,value)
%
%  plot level surface where function is equal to value
%
%  x,y,z:    vectors specifying grid for x,y,z values
%  U    :    3 dimensional array of function values at 3D grid specified by x,y,z
%            
%  use meshgrid to compute U:
%  [X,Y,Z] = meshgrid(x,y,z);  u =  expression in X,Y,Z using .* ./ .^
%
%isosurf(x,y,z,U,value,color)
%    
%  use color for surface (default is to color by value, using caxis)
%     color:   'r', 'g', 'b', 'c', 'm', 'y', 'k'
%            or 'red','green','blue','cyan','magenta','yellow','black'
%            or [r g b]: e.g. [1 0 0] for red
%
%Example: plot surface where x^2 - y^2 - z^2 = 1  for x,y,z in [-3,3]
%
%  x = -3:.2:3; y = -3:.2:3; z = -3:.2:3;  % specify x,y,z grid
%  [X,Y,Z] = meshgrid(x,y,z);
%  U = X.^2 - Y.^2 - Z.^2;        % compute function at grid points
%  isosurf(x,y,z,U,1)             % plot level surface where U = 1

newplot

if nargin==5
  c = 'interp';
  p = patch(isosurface(x,y,z,u,v,v*ones(size(u))));
else
  p = patch(isosurface(x,y,z,u,v));
end
isonormals(x,y,z,u,p);
set(p, 'FaceColor', c, 'EdgeColor', 'none');

daspect([1 1 1]);
view(3);
camlight; lighting phong;
