Let $$b^2h^2m^4-2b^2hkm^3+ (a^2h^2+b^2k^2-a^4e^4)m^2-2a^2hkm+a^2k^2=0$$ with $e^4=\left(-\frac{b^2-a^2}{a^2}\right)^2$
Let $m_1,m_2,m_3,m_4$ be the four solutions. Let $m_1,m_2$ be the slopes of mutually perpendicular normals to the ellipse. So the sums and products of roots of this quartic equations are as follows:-
$$\begin{align} m_1+m_2+m_3+m_4&=\frac{2k}{h}\\ m_1m_2+m_1m_3+m_1m_4+m_2m_3+m_2m_4+m_3m_4&=\frac{a^2h^2+b^2k^2-a^4e^4}{b^2h^2}\\ m_1m_2m_3+m_1m_2m_4+m_1m_3m_4+m_2m_3m_4&=\frac{2a^2k}{b^2h}\\ m_1m_2m_3m_4&=\frac{a^2k^2}{b^2h^2}\end{align}$$
and $$m_1m_2=-1$$
Now how to eliminate $m_1,m_2,m_3,m_4$ from these five equations, so that we get the equation in terms of $a,b,h,k$?
How to derive $(a^2+b^2)(h^2+k^2)(a^2k^2+b^2h^2)^2=(a^2-b^2)^2(a^2k^2-b^2h^2)^2\tag{1}$ from these above five equations?
Now, I got another equations also after doing some algebra.
$(b^2h^2-a^2k^2)m_1^2-(2a^2k+2b^2k)m_1+a^2k^2-b^2h^2=0 \tag{2}$
The roots of equation(3) are here
$m_2=-\frac{1}{m_1},m_3=\frac{a*k}{b*h},m_4=-\frac{a*k}{b*h}$
But these roots are not satifying the following equation. $$m_1m_2m_3+m_1m_2m_4+m_1m_3m_4+m_2m_3m_4=\frac{2a^2k}{b^2h}$$
The above equation can be written as $$-(m_3+m_4)-\frac{2a^2k^2}{b^2h^2}*(m_1+m_2)=\frac{2a^2k}{b^2h^2}\tag{3}$$
Since $m_3+m_4=0, m_1+m_2=\frac{2k}{h}$
Equation (3) becomes $-\frac{2a^2k^3}{b^2h^3}=\frac{2a^2k}{b^2h}$
Does it mean$-h=k$ and/or vice versa? Would any member explain me this discrepancy?
I also got $(a^2-b^2)^2=(a^2+b^2)(h^2+k^2)$ which will be helpful in arriving at equation(1).
How to plot this equation in www.wolframalpha.com?
Solving this problem is possible by hand, but takes a long time. The exact process of solving a quartic equation is outlined here.
Once you solve the quartic equation, you will get the values of $m_1$, $m_2$, $m_3$, and $m_4$ in terms of $a,b,h,$ and $k$. Then, you substitute those values into these equations, giving you what you wanted.