One way is the extended Euclidean algorithm. This can be done in analogy to Gaussian elimination or triangularization, except that, since the coefficient ring is not a field, one has to use the division / Euclidean algorithm to to iteratively descrease the coefficients till zero. Augmenting (appending) an identity matrix will accumulate the effect of the elementary operations, hence yielding the sought Bezout equation. Below is an example from one of my old posts. It computes the Bezout representation for $\rm\:gcd(80,62) = 2\ $ viz. $\ 7\cdot 80\: -\: 9\cdot 62\ =\ 2\:.\:$
For example, to solve m x + n y = gcd(m,n) one begins with two rows [m 1 0], [n 0 1], representing the two equations m = 1m + 0n, n = 0m + 1n. Then one executes the Euclidean algorithm on the numbers in the first column, doing the same operations in parallel on the other columns, Here is an example: d = x(80) + y(62) proceeds as: in equation form | in row form ---------------------+------------ 80 = 1(80) + 0(62) | 80 1 0 62 = 0(80) + 1(62) | 62 0 1 row1 - row2 -> 18 = 1(80) - 1(62) | 18 1 -1 row2 - 3 row3 -> 8 = -3(80) + 4(62) | 8 -3 4 row3 - 2 row4 -> 2 = 7(80) - 9(62) | 2 7 -9 row4 - 4 row5 -> 0 = -31(80) -40(62) | 0 -31 40 Above the row operations are those resulting from applying the Euclidean algorithm to the numbers in the first column, row1 row2 row3 row4 row5 namely: 80, 62, 18, 8, 2 = Euclidean remainder sequence | | for example 62-3(18) = 8, the 2nd step in Euclidean algorithm becomes: row2 -3 row3 = row4 on the identity-augmented matrix.
In effect we have row-reduced the first two rows to the last two. The matrix effecting the reduction is in the bottom right corner. It starts as 1, and is multiplied by each elementary row operation, hence it accumulates the product of all the row operations, namely:
$ \left[ \begin{array}{ccc} 7 & -9\\ -31 & 40\end{array}\right ] \left[ \begin{array}{ccc} 80 & 1 & 0\\ 62 & 0 & 1\end{array}\right ] \ =\ \left[ \begin{array}{ccc} 2 & 7 & -9\\ 0 & -31 & 40\end{array}\right ]\qquad\qquad\ $
Notice row 1 is the particular solution 2 = 7(80) - 9(62) Notice row 2 is the homogeneous solution 0 = -31(80) + 40(62), so the general solution is any linear combination of the two: n row1 + m row2 -> 2n = (7n-31m) 80 + (40m-9n) 62 The same row/column reduction techniques tackle arbitrary systems of linear Diophantine equations. Such techniques generalize easily to similar coefficient rings possessing a Euclidean algorithm, e.g. polynomial rings F[x] over a field, Gaussian integers Z[i]. There are many analogous interesting methods, e.g. search on keywords: Hermite / Smith normal form, invariant factors, lattice basis reduction, continued fractions, Farey fractions / mediants, Stern-Brocot tree / diatomic sequence.