function abc
  x=zeros(21,21);
  y=zeros(21,21);
  z=zeros(21,21);
  x=-1:.1:1;
  y=-1:.1:1;
  for i=1:21
    for j=1:21
      x(i,j)=(j-11)/10.0;
      y(i,j)=(i-11)/10.0;
      z(i,j)=a1(x(i,j),y(i,j));
    end
  end
  figure(1);
  surf(x,y,z)
  title('a1(x,y)');
  for i=1:21
    for j=1:21
      z(i,j)=b1(x(i,j),y(i,j));
    end
  end
  figure(2)
  surf(x,y,z)
  title('b1(x,y)');
  for i=1:21
    for j=1:21
      z(i,j)=c1(x(i,j),y(i,j));
    end
  end
  figure(3)
  surf(x,y,z)
  title('c1(x,y)');
  for i=1:21
    for j=1:21
      z(i,j)=b1(x(i,j),y(i,j))*b1(x(i,j),y(i,j))-a1(x(i,j),y(i,j))*c1(x(i,j),y(i,j));
    end
  end
  figure(4)
  surf(x,y,z)
  title('b1(x,y)^2 - a1(x,y) * c1(x,y)');
  hold on;
  for i=1:21
    for j=1:21
      z(i,j)=0.0;
    end
  end
  surf(x,y,z)
return
  function aa=a1(xx,yy)
    aa=exp(xx/2.0).*exp(yy)/2.0;
  return
  end
  function bb=b1(xx,yy)
    bb=0.7/((xx.*xx.*yy.*yy)+0.5);
  return
  end
  function cc=c1(xx,yy)
    cc=(4.0-exp(xx)-exp(yy/2.0)).*2;
  return
  end
end