Move your /. inside your NSolve, provide a value for y0, add an ==0 or something similar to the end of your second equation, fix several [] that possibly should be (), discard the unused & and see if there are any more errors. Perhaps something like this:

Yes, but do not converge.

That rapidly provides some points very close to solutions of your original system. Perhaps this will be good enough for you. If not then perhaps you can use those as starting locations to then look for more precise solutions, but if you inspect all the results you can see some values of i find very large minima.

Bill,

I notice you square each equation, is there any reason to do that?

Bill,

I suppose  that, as I saw your code. 
In fact, these two expressions, are originally form a gradient (2D )apllied over an energy equation.
The first result expression is concern to the i-componente and the second one to the j-component.

@Cicco, 

Since these equations are rather complicated, I would use some graphical method to handle them first before I do the computation: 

I can quickly check the results, typically the distribution of the root from the following animation: 

Lets take a closer look at the small region around intersection 

Now I have some intuition about where the root might be (actually root loci ). If we are just looking for one possible result, it should be straightforward. For accuracy and efficiency, I would recommend using the optimization method. Because we are looking for the instances that both expressions equal zero, it is equivalent to find the minimum value of the square sum of  the pair of expressions, which ideally should be zero (aka least square = 0). Thus the corresponding stationary point would be one option of the roots. Particularly, FindMinimum function is the one that I would use.  

To have the optimization function return a precise result, it is usually important to set a good initial point. The animation tells us that the diagnal line is not a bad choice. Also as the parameter i moves from 0 to ten, one intersection tends to move from 0 to -10. Thus I can track the root with this piece of code: (Chop clears the residues) 

The animation of contours of the cost functions is illustrated below (the white dot is the root/ the darker the region is, the smaller the cost function/square sum is )

Here offset means how much the least squre sum is away from zero. I also heuristically adjust the initial guess to  {-i+0.1, -i+0.1} so I can have a better approximation. If you just want the numeric result, you simply remove the Manipulate wrapper and use a table function to compute "res" over a list.  

Copiable code:

Hi Yang,
After studying your graphic, I did the following code that works, or at least output the list of pairs (x,y):

$HistoryLength=0;
ClearSystemCache[];
ClearAll["Global`*"];
SetDirectory["/home/marcelo/Documentos/Tomlinson/2D"];
Module[{x,y,i,j, FmolaX,FmolaY,c={0,0},a=2.46,k=1.5625,V=0.5},

  
Emin[go_,ho_]:=Module[{g=go,h=ho},
{x,y}/.Take[Flatten[N[FindMinimum[-V *(2 Cos[(2 Pi x)/a]* Cos[(2 Pi y)/(a Sqrt[3])]+Cos[(4 Pi y)/(a Sqrt[3])])+ (k/2)*(x-g)^2+(k/2)*(y-h)^2,{{x,First},{y,Last}}
,Method->"PrincipalAxis"],7]],-2]];

                                       Do[ 
                                            If[j==0,c={0,i},c={First,Last}];
                                            c=Emin[j,i]  ;                            
                                             FmolaX=-16*k*(First-(j))*(10^-10);
                                             FmolaY=-16*k*(Last-(i))*(10^-10);
                                             {FmolaX,FmolaY}>>>"DataFxFyTomlinsonV15.txt";
                                              c>>>"DataMinimTomlinsonV15.txt",
                                            
                                         (*Print["the minima of the function related to step  Y",i," X",j,":",c],*)  {i,0.0,10.0,0.01},
                                                  {j,0.0,10.0,0.01}]
                                       ];


I explain the moto:

the equation :

-V *(2 Cos[(2 Pi x)/a]* Cos[(2 Pi y)/(a Sqrt[3])] +
    Cos[(4 Pi y)/(a Sqrt[3])]) + (k/2)*(x - g)^2 + (k/2)*(y - h)^2


has to be minimized, which results  {x,y}, I need them as input for the next step of my "Do", they are insert as g0 , h0, respectively: Emin[g0,h0] function.
There is a detail to be respect: as soon as the loop reach (i, "10")the next step (i,"0") do not should get the previous x,y, so I use the "If "command to prevent this.

See that I am runnig across i: 0 to 10, thgrouh J: 0 to 10, at each step ofI I run J :0 to 10, ok?

Important:   "i" corespond to h0 and "J " corresppond to g0     ( x0 is g0 and h0 is  y0, for each step -(x-x0)² and (y-yo)²)

But when I graphic the results (as example the DataFxFyTomlinsonV15.txt ,results) after doing some mining data

Unprotect;
ListDensityPlot[{elemFx}, Mesh -> None, InterpolationOrder -> 2,
ColorFunction -> GrayLevel]
Unprotect;


The picture is not correct to the pic that I have to reproduce....

The core is I need always the very next point (x,y) minimum, for the following step minimizing.

Are my codes correct?

Hi Moderation,
I did the above instructions, (I was using, Open Suse 13.1, kernel version 13.10. , and firefox late version), but every time I published my reply, my text and codes became a terrible mess, I tried several times, and finally I had to use crt C crt V, to post my comments.
Probably a bug?