A trick I find often helps with ellipses is to think of them as stretched/rotated copies of the unit circle.
By reversing the transformation, we can map the problem back onto the unit circle, where the math is usually easier, then run the answer back through the transformation to answer our original question about the ellipse.
To roughly outline the steps:
Find the transformation matrix that maps the unit circle to your ellipse, then compute its inverse
Apply its inverse to vectors along the horizontal and vertical axes (1, 0) & (0, 1) - these give you the directions of the edges of the blue bounding box, after it's been squashed to a diamond.
Our extreme points in our original coordinate system map to where these slanted lines are tangent to the circle.
Take the perpendiculars of these transformed directions (green lines in the diagram above) and normalize them to unit length to snap them onto the unit circle. These are the points of tangency.
Apply your transformation matrix to these tangent points to map them back to the extremes in our original coordinate system.
Our extreme points in our original coordinate system map to where these slanted lines are tangent to the circle.
Take the perpendiculars of these transformed directions (green lines in the diagram above) and normalize them to unit length to snap them onto the unit circle. These are the points of tangency.
Apply your transformation matrix to these tangent points to map them back to the extremes in our original coordinate system.
