It baffles me a bit that people rush to iterative methods, even though in this case there are very straightforward direct methods. A conic section will implicitly be described by a degree $2$ algebraic equation of the form
$$Ax^2+Bxy+Cy^2+Dx+Ey+F=0$$
so if you are given a few points $(x_1,y_1),(x_2,y_2),...$ for each point you get an equation
$$Ax_1^2+Bx_1y_1+Cy_1^2+Dx_1+Ey_1+F=0$$
$$Ax_2^2+Bx_2y_2+Cy_2^2+Dx_2+Ey_2+F=0$$
they are $\textit{linear}$ in the coefficients $A,B,C,...$ of the conic section. So given enough non-degenerate points (technical term "in general poisiton") this homogenous system will have a unique solution. 'Enough' here means five since we will only know $A,B,C,...$ up to a scale factor. If you have more than five points, the system of equations will become overdetermined.
As a joke i once implemented this as a java applet - and it's still up. Maybe this will convince you that this method runs at 'interactive speed' (for five points at least...).
Of course the hard bit then becomes to work with the conic section in the form of the coefficients $A,B,C,...$ - so maybe you have to convert it to the usual parametric form. For this i refer you to for example the ellipse page of math world (equation 15-22). But wikipedia also has more on this.
A sure way to determine whether a point lies on the ellipse is to substitute the point's x- and y-coordinates into the equation and see whether the equation is exactly satisfied. (Your method is also valid, but requires more work.)
An example of how to use MATLAB to plot a curve and points is at https://www.mathworks.com/help/curvefit/fit.html. You will need to solve your equation for y. Be careful, though, because a point may appear to be on the ellipse but actually not be on it.
Best Answer
Given the real parameters $a$, $b$, $c$, $d$, $e$, $f$, the Cartesian equation:
$$ a\,x^2 + b\,x\,y + c\,y^2 + d\,x + e\,y + f = 0 $$
represents a generic conic in the $\left \langle x,\,y \right \rangle$ plane.
In particular, defined the cubic, quadratic and linear invariants:
$$ i \equiv -\frac{1}{2}\det \begin{bmatrix} 2\,a & b & d \\ b & 2\,c & e \\ d & e & 2\,f \\ \end{bmatrix}, \quad \quad j \equiv -\det \begin{bmatrix} 2\,a & b \\ b & 2\,c \\ \end{bmatrix}, \quad \quad k \equiv \text{tr} \begin{bmatrix} 2\,a & b \\ b & 2\,c \\ \end{bmatrix} $$ if: $$ j < 0 \quad \quad \text{and} \quad \quad i\,k > 0 $$ this conic turns out to be a real non-degenerate ellipse.
In this case, a natural parameterization is as follows:
$$ \begin{cases} x = A + B\,\cos u + C\,\sin u \\ y = D + E\,\cos u + F\,\sin u \\ \end{cases} \quad \quad \text{with} \; u \in [0,\,2\,\pi) $$
where is it:
$$ \begin{aligned} & A \equiv \frac{2\,c\,d - b\,e}{j}\,, \quad \quad B \equiv -\frac{2\,\sqrt{c\,i}}{j}\,, \quad \quad C \equiv 0\,, \\ & D \equiv \frac{2\,a\,e - b\,d}{j}\,, \quad \quad E \equiv \frac{b}{j}\,\sqrt{\frac{i}{c}}\,, \quad \quad F \equiv \sqrt{\frac{-i}{c\,j}}\,. \end{aligned} $$
At this point it's necessary to make it clear what we're interested in.
If we were interested in the coordinates of the vertices of this ellipse, they're determined by minimizing and maximizing the distance function between a generic point of the ellipse and the center of the ellipse, i.e.
if $b = 0$, the two pairs of vertices are identified by the following angles: $$ \begin{aligned} & u_1 = 0\,, \quad \quad u_2 = \pi\,; \\ & u_3 = \frac{\pi}{2}\,, \quad \quad u_4 = \frac{3\,\pi}{2}\,; \end{aligned} $$
if $b \ne 0$, the two pairs of vertices are identified by the following angles: $$ \begin{aligned} & u_1 = \arctan(z^+)\,, \quad \quad u_2 = \arctan(z^+) + \pi\,; \\ & u_3 = \arctan(z^-) + \pi\,, \quad \quad u_4 = \arctan(z^-) + 2\,\pi\,; \end{aligned} $$
where is it:
$$ z^{\pm} \equiv \frac{-\left(B^2-C^2+E^2-F^2\right) \pm \sqrt{\left(B^2-C^2+E^2-F^2\right)^2 + 4\,\left(B\,C + E\,F\right)^2}}{2\left(B\,C + E\,F\right)}\,. $$
Instead, if we were interested in the coordinates of the extreme points of this ellipse, they're determined by minimizing and maximizing respectively the abscissa and the ordinate of a generic point of the ellipse, i.e. $$ u_1 = 0\,, \quad \quad u_2 = \pi\,; \quad \quad u_3 = 2\arctan(z^+)\,, \quad \quad u_4 = 2\arctan(z^-) + 2\,\pi\,; $$
where is it:
$$ z^{\pm} \equiv \frac{-E \pm \sqrt{E^2 + F^2}}{F}\,. $$
I wish you a good study.