Naturally linear Diophantine equations
Pyth, 92 bytes

I!%vzhK%2u?sm,ed-hd*ed/F&lt;G2cG2@G1G+~Q,hQ_eQj9 2)J*L/vzhKtKeoSNm-VJ/RhK_*LdQsm+LdtM3/V*LhK_JQ


It's quite a monster. 

Try it online: Demonstration. The input format is c\n[a,b] and the output format is [x,y]. 

In the case that no integer solution exists, I'll print nothing, and in the case that no natural integer solution exists, I'll simply print a random integer solution. 

Explanation (Rough Overview)


At first I'll find an integer solution to the equation ax + by = gcd(a,b) by using the Extended Euclidean algorithm. 
Then I'll modify the solution (my multiplying a and b with c/gcd(a,b)) to get an integer solution of ax + by = c. This works, if c/gcd(a,b) is an integer. Otherwise there doesn't exist a solution. 
All the other integer solutions have the form a(x+n*b/d) + b(y-n*a/d) = c  with d = gcd(a,b) for integer n. Using the two inequalities x+n*b/d &gt;= 0 and y-n*a/d &gt;= 0 I can determine 6 possible values for n. I'll try all 6 of them and print the solution with the highest lowest coefficient. 


Explanation (Detailed)

The first step is to find an integer solution to the equation ax' + by' = gcd(a,b). This can be done by using the extended Euclidean algorithm. You can get an idea on how it works at Wikipedia. The only difference is, that instead of using 3 columns (r_i s_i t_i) I'll use 6 columns (r_i-1 r_i s_i-1 s_i t_i-1 t_i). This way I don't have to keep the last two rows in memory, only the last one. 

K%2u?sm,ed-hd*ed/F&lt;G2cG2@G1G+~Q,hQ_eQj9 2)   implicit: Q = [a,b] (from input)
                                     j9 2    convert 9 to base 2: [1,0,0,1]
                            + Q              add to Q =&gt; [a,b,1,0,0,1]
                                             this is the initial row
   u                                     )   start with G = ^ and update G repeatedly
                                             by the following expression, until
                                             the value of G doesn't change anymore
    ?                   @G1                    if G[1] != 0:
                     cG2                         split G into parts of 2
      m                                          map the parts d to:
       ,                                           the pair 
        ed                                           d[1]
          -hd*ed/F&lt;G2                                d[0]-d[1]*G[0]/G[1]
     s                                           unfold
                                               else:
                           G                     G (don't change it, stop criterion for u)
 %2                                          take every second element
                                             we get the list [gcd(a,b),x',y']
K                                            store this list in K
                             ~Q,hQ_eQ        afterwards change Q to [Q[0],-Q[1]] = [a,-b]
                                             This will be important for the other parts. 


Now I want to find a solution to ax + by = c. This is possible only, when c mod gcd(a,b) == 0. If this equation is satisfied, I simply multiplying x',y' with c/gcd(a,b). 

I!%vzhK...J*L/vzhKtK   implicit: z = c in string format (from input)
  %vzhK                evaluated(z) mod K[0] (=gcd(a,b))
I!                     if not ^ than: 
             /vzhK        c/K[0]
           *L     tK      multipy ^ to each element in K[1:] (=[x',y'])
          J               and store the result in J, this is now [x,y]


We have an integer solution for ax + by = c. Notice, that x, y or both may be negative. So our goal is to transform these to non-negative. 

The nice thing about Diophantine equations is, that we can describe all solution using only one initial solution. If (x,y) is a solution, that all other solutions are of the form (x-n*b/gcd(a,b),y+n*a/gcd(a,b)) for n integer. 

Therefore we want to find a n, where x-n*b/gcd(a,b) &gt;= 0 and y+n*a/gcd(a,b &gt;= 0. After some transformation we end up with the two inequalities n &gt;= -x*gcd(a,b)/b and n &gt;= y*gcd(a,b)/a. Notice that the inequality symbol might look in the other direction due the division with a potential negative a or b. I don't care that much about it, I simply say that one number of -x*gcd(a,b)/b - 1, -x*gcd(a,b)/b, -x*gcd(a,b)/b + 1 definitly satisfies inequality 1, and one number of y*gcd(a,b)/a - 1, y*gcd(a,b)/a, y*gcd(a,b)/a + 1 satisfies inequality 2. It there is a n, that satisfies both inequalities, one of the 6 numbers also does. 

Then I calculate the new solutions (x-n*b/gcd(a,b),y+n*a/gcd(a,b)) for all 6 possible values of n. And I print the solution with the highest lowest value. 

eoSNm-VJ/RhK_*LdQsm+LdtM3/V*LhK_JQ
                               _J    reverse J =&gt; [y,x]
                           *LhK      multiply each value with K[0] =&gt; [y*gcd,x*gcd]
                         /V      Q   vectorized division =&gt; [y*gcd/a,-x*gcd/b]
                  m                  map each d of ^ to:
                      tM3              [-1,0,1]
                   +Ld                 add d to each ^
                 s                   unfold
                                     these are the possible values for n
    m                                map each d (actually n) of ^ to:
             *LdQ                      multiply d to Q =&gt; [a*n,-b*n]
            _                          reverse =&gt; [-b*n,a*n]
        /RhK                           divide by K[0] =&gt; [-b*n/gcd,a*n/gcd]
     -VJ                               vectorized subtraction with J
                                       =&gt; [x+b*n/gcd,y-a*n/gcd]
 oSN                                 order the solutions by their sorted order
e                                    print the last one


The sort by their sorted order thing works the following way. I'm using the example 2x + 3y = 11

I sort each of the 6 solutions (this are called keys), and sort the original solutions by their keys:

solutions: [1, 3], [4, 1], [7, -1], [-5, 7], [-2, 5], [1, 3]
keys:      [1, 3], [1, 4], [-1, 7], [-5, 7], [-2, 5], [1, 3]
sort by key:
solutions: [-5, 7], [-2, 5], [7, -1], [1, 3], [1, 3], [4, 1]
keys:      [-5, 7], [-2, 5], [-1, 7], [1, 3], [1, 3], [1, 4]


This sorts a complete non-negative solution to the end (if there is any).