I want to find a set of integer solutions of Diophantine equation: $ax + by = c$, and apparently $\gcd(a,b)|c$. Then by what formula can I use to find $x$ and $y$ ? 

I tried to play around with it:
$x = (c - by)/a$, hence $a|(c - by)$.   

$a$, $c$ and $b$ are known. So to obtain integer solution for $a$, then $c - by = ak$, and I lost from here, because $y = (c - ak)/b$. I kept repeating this routine and could not find a way to get rid of it? Any hint?

Thanks,
Chan

The diophantine equation $ax+by = c$ has solutions if and only if $\gcd(a,b)|c$. If so, it has infinitely many solutions, and any one solution can be used to generate all the other ones.

To see this, note that the greatest common divisor of $a$ and $b$ divides both $ax$ and $by$, hence divides $c$ if there is a solution. This gives the necessity of the condition (which you have backwards). (fixed in edits)

The converse is actually a constructive proof, that you can find in pretty much every elementary number theory course or book, and which is essentially the same as yunone's answer above (but without dividing through first).

From the Extended Euclidean Algorithm, given any integers $a$ and $b$ you can find integers $s$ and $t$ such that $as+bt = \gcd(a,b)$; the numbers $s$ and $t$ are not unique, but you only need one pair. Once you find $s$ and $t$, since we are assuming that $\gcd(a,b)$ divides $c$, there exists an integer $k$ such that $\gcd(a,b)k = c$. Multiplying $as+bt=\gcd(a,b)$ through by $k$ you get
$$a(sk) + b(tk) = \gcd(a,b)k = c.$$

So this gives one solution, with $x=sk$ and $y=tk$. 

Now suppose that $ax_1 + by_1 = c$ is a solution, and $ax+by=c$ is some other solution. Taking the difference between the two, we get
$$a(x_1-x) + b(y_1-y) = 0.$$
Therefore, $a(x_1-x) = b(y-y_1)$. That means that $a$ divides $b(y-y_1)$, and therefore $\frac{a}{\gcd(a,b)}$ divides $y-y_1$. Therefore, $y = y_1 + r\frac{a}{\gcd(a,b)}$ for some integer $r$. Substituting into the equation $a(x_1-x) = b(y-y_1)$ gives
$$a(x_1 - x) = rb\left(\frac{a}{\gcd(a,b)}\right)$$
which yields
$$\gcd(a,b)a(x_1-x) = rba$$
or $x = x_1 - r\frac{b}{\gcd(a,b)}$. 

Thus, if $ax_1+by_1 = c$ is any solution, then all solutions are of the form
$$x = x_1 - r\frac{b}{\gcd(a,b)},\qquad y = y_1 + r\frac{a}{\gcd(a,b)}$$
exactly as yunone said.


To give you an example of this in action, suppose we want to find all integer solutions to $$258x + 147y = 369.$$

First, we use the Euclidean Algorithm to find $\gcd(147,258)$; the parenthetical equation on the far right is how we will use this equality after we are done with the computation.
\begin{align*}
258 &= 147(1) + 111 &\quad&\mbox{(equivalently, $111=258 - 147$)}\\
147 &= 111(1) + 36&&\mbox{(equivalently, $36 = 147 - 111$)}\\
111 &= 36(3) + 3&&\mbox{(equivalently, $3 = 111-3(36)$)}\\
36 &= 3(12).
\end{align*}
So $\gcd(147,258)=3$. Since $3|369$, the equation has integral solutions.

Then we find a way of writing $3$ as a linear combination of $147$ and $258$, using the Euclidean algorithm computation above, and the equalities on the far right. We have:
\begin{align*}
3 &= 111 - 3(36)\\
&= 111 - 3(147 - 111) = 4(111) - 3(147)\\
&= 4(258 - 147) - 3(147)\\
&= 4(258) -7(147).
\end{align*}
Then, we take $258(4) + 147(-7)=3$, and multiply through by $123$; why $123$? Because $3\times 123 = 369$. We get:
$$258(492) + 147(-861) = 369.$$
So one solution is $x=492$ and $y=-861$. All other solutions will have the form
\begin{align*}
x &= 492 - \frac{147r}{3} = 492 - 49r,\\
y &= -861 + \frac{258r}{3} =86r - 861, &\qquad&r\in\mathbb{Z}.
\end{align*}
You can reduce those constants by making a simple change of variable. For example, if we let $r=t+10$, then 
\begin{align*}
x &= 492 - 49(t+10) = 2 - 49t,\\
y &= 86(t+10) - 861 = 86t - 1,&\qquad&t\in\mathbb{Z}.
\end{align*}
