Solution 1:

We claim that every integer $z \geq 49x^2y^2$ is in $S_{x,y}$. Since $\gcd(x,y)=1$, there exist $u,v\in\mathbb{Z}$ such that $ux-vy=1$. We may assume that $0< u\leq y$ and $0\leq v < x$. Now, note that the interval $I:=\left(\frac{zv}{x},\frac{zu}{y}\right)$ is of measure $$\frac{zu}{y}-\frac{zv}{x}=\frac{z}{xy}(ux-vy)=\frac{z}{xy}\geq 7\sqrt{z}>1+6\sqrt{z}\,.$$ Consequently, $I$ contains at least $\lceil6\sqrt{z}\rceil$ consecutive integers, whence $I$ contains an integer $p$ such that $\gcd(p,z)=1$, by the lemma below. Note that $\frac{zv}{x}<p<\frac{zu}{y}$.

To solve $n_1x+n_2y=z$ for $n_1,n_2\in\mathbb{N}$, we consider solutions $n_1=zu-p y>0$ and $n_2=px-zv>0$. We claim that $\gcd\left(n_1,n_2\right)=1$. First, from the assumption that $\gcd(p,z)=1$, there are $\mu,\nu\in\mathbb{Z}$ such that $p\mu+z\nu=1$. Now, take $a:=\mu v+\nu x$ and $b:=\mu u+\nu y$. We see that $$\begin{align} an_1+bn_2 &=(\mu v+\nu x)\big(zu-py\big)+(\mu u+\nu y)\big(px-zv\big) \\ &=z\big(u(\mu v+\nu x)-v(\mu u+\nu y)\big)+p\big(x(\mu u+\nu y)-y(\mu v+\nu x)\big) \\ &=\nu(ux-vy)z+\mu(ux-vy)p=\mu p +\nu z=1\,, \end{align}$$ which implies that $\gcd\left(n_1,n_2\right)=1$. Hence, $z\in S_{x,y}$.

Lemma: For every positive integer $N$, a set $L$ composed by consecutive integers such that $|L|\geq\lceil6\sqrt{N}\rceil$ must contain an integer $t$ which is coprime to $N$.

(The lemma can be generalized. For every $\epsilon>0$, there exists a constant $\kappa_\epsilon>0$ for which the following holds: for any $N \in \mathbb{N}$, $d\in\mathbb{N}$, and $r\in\mathbb{Z}$ such that $\gcd(N,d)=1$, a set $L$ composed by consecutive elements of the arithmetic progression $\big\{nd+r\big\}_{n\in\mathbb{Z}}$ with $|L| \geq \left\lceil \kappa_\epsilon N^\epsilon\right\rceil$ has an element $t\in L$ such that $\gcd(t,N)=1$. If $\epsilon \geq 1$, then trivially, $\kappa_\epsilon$ can be taken to be $1$. A proof for this generalized version with $0<\epsilon<1$ is done similarly to the one given below.)

Proof: Suppose that $N$ has (pairwise distinct) prime factors $p_1,p_2,\ldots,p_r \in\mathbb{N}$. Let $L$ be a set composed by $l$ consecutive integers, where $l\geq \lceil 6\sqrt{N}\rceil$.

Let $[r]:=\{1,2,\ldots,r\}$. For each subset $T$ of $[r]$, define $f_L(T)$ to be the number of elements of $L$ divisible by $p_j$ for all $j\in T$. Obviously, $\left\lfloor \frac{l}{\prod_{\lambda \in T} p_\lambda}\right\rfloor \leq f_L(T) \leq \left\lceil \frac{l}{\prod_{\lambda \in T} p_\lambda}\right\rceil$. Hence, the number of elements of $L$ which are coprime to $N$ is given by $$K_L:=\sum_{T\subseteq[r]}\,(-1)^{|T|}f_L(T)\,,$$ where we have used the Principle of Inclusion and Exclusion.

If $T\subseteq[r]$ is such that $|T|$ is even, $$f_L(T)> \frac{l}{\prod_{\lambda \in T} p_\lambda}-1=\frac{l}{\prod_{\lambda \in T} p_\lambda}-(-1)^{|T|}\,,$$ if $T\subseteq[r]$ is such that $|T|$ is odd, $$f_L(T)<\frac{l}{\prod_{\lambda \in T} p_\lambda}+1=\frac{l}{\prod_{\lambda \in T} p_\lambda}-(-1)^{|T|}\,.$$ This means $$ \begin{align} K_L &>\sum_{T\subseteq[r]}\,\left(\frac{(-1)^{|T|}l}{\prod_{\lambda \in T} p_\lambda}-1\right)=l\,\sum_{T\subseteq[r]}\frac{(-1)^{|T|}}{\prod_{\lambda\in T} p_\lambda}-\sum_{T\subseteq[r]}1 \\ &=l\,\prod_{j=1}^r\,\left(1-\frac{1}{p_j}\right)-2^r \geq 6\sqrt{N}\,\prod_{j=1}^r\,\left(1-\frac{1}{p_j}\right)-2^r \\ & \geq 6\sqrt{\prod_{j=1}^r\,p_j}\,\prod_{j=1}^r\,\left(1-\frac{1}{p_j}\right)-2^r=6\,\prod_{j=1}^r\,\left(\frac{p_j-1}{\sqrt{p_j}}\right)-2^r\,. \end{align}$$

If we order the positive primes as $q_1<q_2<q_3<\ldots$, then $$6\,\prod_{j=1}^r\,\left(\frac{p_j-1}{\sqrt{p_j}}\right) \geq 6\,\prod_{j=1}^r\,\left(\frac{q_j-1}{\sqrt{q_j}}\right)\,.$$ Observe that, for an integer $j \geq 4$, $\frac{q_j-1}{\sqrt{q_j}}>2$. By induction on $r$ (where the base cases $r=1,2,3$ must be manually checked), we conclude that $6\,\prod_{j=1}^r\,\left(\frac{q_j-1}{\sqrt{q_j}}\right)>2^r$ for every $r \in \mathbb{N}$. Thus, $$K_L > 6\,\prod_{j=1}^r\,\left(\frac{p_j-1}{\sqrt{p_j}}\right)-2^r \geq 6\,\prod_{j=1}^r\,\left(\frac{q_j-1}{\sqrt{q_j}}\right)-2^r > 0\,.$$ Ergo, $L$ contains an element $t$ such that $\gcd(t,N)=1$.

P.S. I believe that the largest integer $m(x,y)$ that is not in $S_{x,y}$ satisfies $m(x,y) \in \Theta(xy)$. What I know is that $m(x,y)\in\Omega(xy)$ and $m(x,y) \in O\left((xy)^{1+\varepsilon}\right)$ for all $\varepsilon>0$. However, it seems like an extremely difficult problem to find an explicit or an asymptotic formula for $m(x,y)$.

EDIT: After some experiments, it seems like $m(x,y)\in\Theta\big((xy)\,\ln(xy)\big)$. A plot to illustrate this conjecture is given below. Here, $x,y\in\{1,2,\ldots,50\}$. The scattered plot in blue is the actual $xy$ versus $m(x,y)$. The red curve is the estimate $m(x,y)=(xy)\,\ln(xy)$. If you have an efficient way to write a code to determine the relationship between $xy$ and $m(x,y)$, please share it with us. My code takes quite long to run.enter image description here