poe := proc(a,b)
       local bAB, bAC, M, Mx, My, triang, circ, rad;
       triang := plot([0,0,1,0,a,b,0,0], style=line, color = green);
       bAB := bisector(0,0,1,0);
       bAC := bisector(0,0,a,b);
       M := solve({bAB, bAC}, {x, y});
       Mx := subs(M, x);
       My := subs(M, y);
       rad := sqrt(Mx^2 + My^2);
       circ := plot([rad*cos(t)+Mx, rad*sin(t)+My, t=0..2*Pi],
                    color = red);
       plot({triang,circ}, axes=framed);
end;

bisector := proc(a,b,c,d)
       (d-b)*(y-(b+d)/2)+(c-a)*(x-(a+c)/2);
end;