I believe the answer should be

```
1600*(9*x/40 + y + 37/40)**2/1681 + 1600*(x - 9*y/40 - 91/10)**2/15129 - 1
```

You can confirm by plotting the equations (in implicit form) at desmos.com/calculator. I used the equation I derived and posted at SO (referenced in issue #2815)

```
>>> ge
-1 + (-s*(x - xc) + y - yc)**2/(m**2*(s**2 + 1)) + (s*(y - yc) + x - xc)**2/(M**2*(s**2 + 1))
>>> e=Ellipse((3,4),1,3)
>>> line = Line(Point(0,-4), Point(5,0))
>>> ge.subs(zip((xc,yc), e.center.reflect(line).args)).subs(s,
... Line(e.center,slope=0 if (e.hradius==e.major) else oo).reflect(line).slope
... ).subs(m,e.minor).subs(M,e.major)
1600*(9*x/40 + y + 37/40)**2/1681 + 1600*(x - 9*y/40 - 91/10)**2/15129 - 1
```