This paper describes a non-linear structure-preserving matrix method for the computation of the coefficients of an approximate greatest common divisor (AGCD) of degree t of two Bernstein polynomials f (y) and g(y). This method is applied to a modified form S t (f, g)Q t of the tth subresultant matrix S t (f, g) of the Sylvester resultant matrix S(f, g) of f (y) and g(y), where Q t is a diagonal matrix of combinatorial terms. This modified subresultant matrix has significant computational advantages with respect to the standard subresultant matrix S t (f, g), and it yields better results for AGCD computations. It is shown that f (y) and g(y) must be processed by three operations before S t (f, g)Q t is formed, and the consequence of these operations is the introduction of two parameters, α and θ, such that the entries of S t (f, g)Q t are non-linear functions of α, θ and the coefficients of f (y) and g(y). The values of α and θ are optimised, and it is shown that these optimal values allow an AGCD that has a small error, and a structured low rank approximation of S(f, g), to be computed.