Here is how I compute a conic without $5\times5$ equations, based on my background in projective geometry.
Finding the matrix
Start with homogeneous coordinates, i.e. you have five points
Now a point $P$ lies on the conic with these five points iff
$$[A,C,E][B,D,E][A,D,P][B,C,P] - [A,D,E][B,C,E][A,C,P][B,D,P] = 0$$
where I use $[\cdot,\cdot,\cdot]$ to denote a determinant. Now you may know that you can write a $3\times3$ matrix as a triple product, e.g.
$$[A,D,P] = \langle A\times D,P\rangle$$
Combine two of these and you have a quadratic form with a rank 1 matrix in the center:
$$[A,D,P][B,C,P] = \langle P,A\times D\rangle\cdot\langle B\times C,P\rangle
= P^T\cdot(A\times D)\cdot(B\times C)^T\cdot P$$
So the original equation boils down to $P^TMP=0$ using the following matrix:
M &=\phantom+ [A,C,E][B,D,E]\cdot(A\times D)\cdot(B\times C)^T \\
&\phantom=- [A,D,E][B,C,E]\cdot(A\times C)\cdot(B\times D)^T
You probably should symmetrize your final result as well, i.e. compute $M+M^T$.
So you have to compute four determinants, four cross products, two outer products, two scalar times matrix products, one matrix subtraction and one matrix addition. But all of the vectors and matrices will be $3\times 3$ only, and you never have to pivot, never have to make any case distinctions. If you work with the homogeneous coordinates using the representatives given above, many numbers in your computations will be equal to $1$, which can be used to further simplify an implementation.
You may notice that the determinants in the left part of each line are just $E$ plugged into the quadratic form you get from the right part of the other line. So if evaluating quadratic forms is any easier for you than computing determinants, go ahead and re-use the matrices you need for the right hand side in any case.
Finding the center
Now you want the center of that beast. The center is the pole of the line at infinity. For that you need the dual matrix, which algebraically is cheapest to compute using the classical adjoint. Multiply that matrix by $(0,0,1)$ and you have the homogeneous coordinates of the center. Divide the first two coordinates by the third to get back to inhomogeneous coordinates.