First multiply through with $a^2b^2$ expand and homogenize the equation:
$$(a^2k^2+b^2h^2-a^2b^2)z^2-2a^2kyz-2b^2hxz+a^2y^2+b^2x^2=0$$
Then the system is
$$\begin{align}X-\lambda (2b^2p-2b^2hr)&=0\\
Y-\lambda (2a^2q-2a^2kr)&=0\\
Z-\lambda (2(a^2k^2+b^2h^2-a^2b^2)r-2a^2kq-2b^2hp)&=0\\
Xp+Yq+Zr&=0\end{align}$$
Solving explicitly is feasible by hand (showing the Sylvester matrix of which the resultant is the determinant):
resultant(X-l*(2*b^2*p-2*b^2*h*r),Y-l*(2*a^2*q-2*a^2*k*r),p);
$$\begin{pmatrix} 2\,a^2\,k
\,l\,r-2\,a^2\,l\,q+Y\\\end{pmatrix}$$
resultant(X-l*(2*b^2*p-2*b^2*h*r),Z-l*(2*(a^2*k^2+b^2*h^2-a^2*b^2)*r-2*a^2*k*q-2*b^2*h*p),p);
$$\begin{pmatrix} -2\,b^2\,
l&2\,b^2\,h\,l\,r+X\\ 2\,b^2\,h\,l&\left(
-2\,a^2\,k^2-2\,b^2\,h^2+2\,a^2\,b^2\right)\,l\,r+2\,a^2\,k\,l\,q+Z
\\ \end{pmatrix}$$
resultant(Y-l*(2*a^2*q-2*a^2*k*r),Z-l*(2*(a^2*k^2+b^2*h^2-a^2*b^2)*r-2*a^2*k*q-2*b^2*h*p),p);
$$\begin{pmatrix} 2\,a^2\,k
\,l\,r-2\,a^2\,l\,q+Y\\\end{pmatrix}$$
resultant(X-l*(2*b^2*p-2*b^2*h*r),X*p+Y*q+Z*r,p);
$$\begin{pmatrix}-2\,b^2\,
l&2\,b^2\,h\,l\,r+X\\ X&Z\,r+Y\,q\\
\end{pmatrix}$$
resultant(Y-l*(2*a^2*q-2*a^2*k*r),X*p+Y*q+Z*r,p);
$$\begin{pmatrix} 2\,a^2\,k
\,l\,r-2\,a^2\,l\,q+Y\\
\end{pmatrix}$$
resultant(Z-l*(2*(a^2*k^2+b^2*h^2-a^2*b^2)*r-2*a^2*k*q-2*b^2*h*p),X*p+Y*q+Z*r,p);
$$\begin{pmatrix} 2\,b^2\,h
\,l&\left(-2\,a^2\,k^2-2\,b^2\,h^2+2\,a^2\,b^2\right)\,
l\,r+2\,a^2\,k\,l\,q+Z\\ X&Z\,r+Y\,q\\
\end{pmatrix}$$
Among these choose two and take resultant w.r.t $q$
resultant((2*X*a^2*k^2+2*X*b^2*h^2+2*Z*b^2*h-2*X*a^2*b^2)*l*r+(2*Y*b^2*h-2*X*a^2*k)*l*q-X*Z,2*a^2*k*l*r-2*a^2*l*q+Y,q);
$$\begin{pmatrix}\left(2\,
Y\,b^2\,h-2\,X\,a^2\,k\right)\,l&\left(2\,X\,a^2\,k^2+2
\,X\,b^2\,h^2+2\,Z\,b^2\,h-2\,X\,a^2\,b^2\right)\,l\,r-X\,Z\\ -2\,a
^2\,l&2\,a^2\,k\,l\,r+Y\\ \end{pmatrix}$$
resultant((-(2*X*b^2*h+2*Z*b^2)*l*r)-2*Y*b^2*l*q-X^2,2*a^2*k*l*r-2*a^2*l*q+Y,q);
$$\begin{pmatrix} -2\,Y\,b^
2\,l&\left(-2\,X\,b^2\,h-2\,Z\,b^2\right)\,l\,r-X^2\cr
-2\,a^2\,l&2\,a^2\,k\,l\,r+Y\\ \end{pmatrix}$$
Finally the resultant w.r.t $r$
resultant(2*l*((2*Y*a^2*b^2*h*k+2*X*a^2*b^2*h^2+2*Z*a^2*b^2*h-2*X*a^4*b^2)*l*r-X*Y*a^2*k+Y^2*b^2*h-X*Z*a^2),-2*l*((b^2*(2*Y*a^2*k+2*Z*a^2)+2*X*a^2*b^2*h)*l*r+Y^2*b^2+X^2*a^2),r);
$$\begin{pmatrix} \left(4\,
Y\,a^2\,b^2\,h\,k+4\,X\,a^2\,b^2\,h^2+4\,Z\,a^2\,b^2\,h-4\,X\,a^4\,b
^2\right)\,l^2&\left(-2\,X\,Y\,a^2\,k+2\,Y^2\,b^2\,h-2
\,X\,Z\,a^2\right)\,l\cr \left(-4\,Y\,a^2\,b^2\,k-4\,X\,a^2\,b^2\,h-
4\,Z\,a^2\,b^2\right)\,l^2&\left(-2\,Y^2\,b^2-2\,X^2\,a
^2\right)\,l\\ \end{pmatrix}$$
factor(%);
-8*X*a^4*b^2*(Y^2*k^2+2*X*Y*h*k+2*Y*Z*k+X^2*h^2+2*X*Z*h-Y^2*b^2-X^2*a^2+Z^2)*l^3
i.e. $Y^2k^2+2XYhk+2YZk+X^2h^2+2XZh-Y^2b^2-X^2a^2+Z^2=0.$
or in one fell swoop using the maxima CAS function eliminate
eliminate([X-l*(2*b^2*p-2*b^2*h*r),Y-l*(2*a^2*q-2*a^2*k*r),Z-l*(2*(a^2*k^2+b^2*h^2-a^2*b^2)*r-2*a^2*k*q-2*b^2*h*p),X*p+Y*q+Z*r],[p,q,r]);
[-16*Y*a^2*b^6*(Y^2*k^2+2*X*Y*h*k+2*Y*Z*k+X^2*h^2+2*X*Z*h-Y^2*b^2-X^2*a^2+Z^2)*l^4]
Setting $Z=1$ we have the dual by the equation in the standard affine in the dual projective plane
$$(h^2-a^2)X^2+2hkXY+(k^2-b^2)Y^2+2hX+2kY+1=0.$$
$$(a^2k^2+b^2h^2-a^2b^2)z^2-2a^2kyz-2b^2hxz+a^2y^2+b^2x^2=0$$
Then the system is
$$\begin{align}X-\lambda (2b^2p-2b^2hr)&=0\ Y-\lambda (2a^2q-2a^2kr)&=0\ Z-\lambda (2(a^2k^2+b^2h^2-a^2b^2)r-2a^2kq-2b^2hp)&=0\ Xp+Yq+Zr&=0\end{align}$$
– Jan-Magnus Økland Mar 15 '21 at 18:37