% find_roots.m  r = roots(p)
function find_roots()
   diary('find_roots_m.out')
   disp('find roots.m running')
   disp('some examples have roots on unit circle')
   disp('x^2+1=0')
   p2 = [1.0 0.0 1.0]
   r2 = roots(p2)
   disp('x^3+1=0')
   p3 = [1.0 0.0 0.0 1.0]
   r3 = roots(p3)
   disp('x^4+1=0')
   p4 = [1.0 0.0 0.0 0.0 1.0]
   r4 = roots(p4)
   disp('x^5+1=0')
   p5 = [1.0 0.0 0.0 0.0 0.0 1.0]
   r5 = roots(p5)
   disp('x^6+1=0')
   p6 = [1.0 0.0 0.0 0.0 0.0 0.0 1.0]
   r6 = roots(p6)
   disp('5 x^2 -6 x +5 =0')
   pq = [5.0 -6.0 5.0]   
   rq = roots(pq)
   disp('x^2 -x + 3  not unit')
   pd = [1.0 -1.0 3.0]
   rd = roots(pd)
   disp('x^2 - 1  real')
   pa = [1.0 0.0 -1.0]
   ra = roots(pa)
   disp('x^2 -3 x + 2  real')
   pr = [1.0 -3.0 2.0]
   rr = roots(pr)
   disp('x^4 -10 x^3 +35 x^2 -50 x + 24  real')
   disp('(x-1)*(x-2)*(x-3)*(x-4)')
   ps = [1.0 -10.0 35.0 -50.0 24.0]
   rs = roots(ps)
   disp('x^2 -50 x + 0.5  big, small')
   pb = [1.0 -50.01 0.5]
   rb = roots(pb)
end
