> restart;

> with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> a:=3:b:=5:r:=2.22:s:=5.1:

> ss:=[a*cos(r),b*sin(r)]:s:=[a*cos(s),b*sin(s)];

> m:=[s,ss];

> for k to 200 do u:=ss-s: n:=[-b*ss[1]/a,-a*ss[2]/b]:v:=u-2*dotprod(u , n)*n/dotprod(n,n):t:=-2*(b^2*v[1]*ss[1]+a^2*ss[2]*v[2])/(b^2*v[1]^2+a^2*v[2]^2):sss:=ss+t*v:m:=[op(m),sss]:s:=ss:ss:=sss:od:

> plot(m,scaling=constrained,axes=none);

>