For $x, y, n\in\mathbb{N}$ find $\lambda, \mu \in \mathbb{N}$ so that $\lambda x + \mu y = n$ We know from the Bézout's identity that for any $x, y \in \mathbb{Z}$ and $d := \gcd(x,y)$ there are $\mu, \lambda \in\mathbb{Z}$ so that 
$$\mu x + \lambda y = d\, .$$
We immediately see that if $d\mid n \in \mathbb{Z}$ we can also find $\mu, \lambda \in\mathbb{Z}$ so that
$$\mu x + \lambda y = n\, .$$
Now this result seems still to be true if we restrict ourselves to natural numbers and $n$ is "large enough":

For $x,y,n \in \mathbb{N}$ and $n$ "large enough" there are $\mu, \lambda \in \mathbb{N}$ so that 
  $$\mu x + \lambda y = n\, .$$



*

*How can this theorem be precisely stated and how would one prove it?

*How can I calculate the $\mu, \lambda$? What does the algorithm look like (maybe its implemented in Sage?)?

 A: Your result is already well stated.
The equation $\mu x + \lambda y = n$ is a linear Diophantine equation, for which you can write out every possible solution (with $\lambda,\mu \in \Bbb{Z}$) once you know one particular solution. Furthermore, finding a particular solution usually isn't hard.
On the other hand, this answer tells you that the equation has either
$$
\left\lfloor \frac{n}{xy} \right\rfloor
\qquad \text{or} \qquad
\left\lfloor \frac{n}{xy} \right\rfloor + 1
$$
solutions with $\lambda,\mu > 0$. In particular, the equation has always at least one solution in positive intgers whenever $n > xy$.

Here's a way to find a positive solution (if your equation admits any): first find a particular solution $(\lambda_0,\mu_0) \in \Bbb{Z}^2$ of your equation. Since every other solution is of the form $\lambda = \lambda_0 + k y/d$, $\mu = \mu_0 - k x/d$ where $d = \gcd(x,y)$ and $k$ is any integer, all that's left to do is solve the system of inequalities
$$
\begin{cases}
\lambda_0 + k \frac{y}{d} > 0 \\
\mu_0 - k \frac{x}{d} > 0
\end{cases}
$$

Note: In SAGE you can easily find all the non-negative solutions by instancing the WeightedIntegerVectors object, while in Mathematica or WolframAlpha you can use the FrobeniusSolve function.
