Let $X=\{x_1,x_2,...,x_n\}$ and $Y=\{y_1,y_2,...,y_n\}$ be sets over a finite field $F$ with $p=char(F)>2$. Assume
$$x_1^k+x_2^k+...+x_n^k=y_1^k+y_2^k+...+y_n^k,\ 1\leq k\leq n$$
I wanna show that $X=Y$.
Let $X=\{x_1,x_2,...,x_n\}$ and $Y=\{y_1,y_2,...,y_n\}$ be sets over a finite field $F$ with $p=char(F)>2$. Assume
$$x_1^k+x_2^k+...+x_n^k=y_1^k+y_2^k+...+y_n^k,\ 1\leq k\leq n$$
I wanna show that $X=Y$.
Use the standard notations $e_k=\sum_{A\subset \{1,\dots,n\}, |A|=k} \prod_{i\in A} x_i$, with the conventions $e_0=1$ and $e_m=0$ for $m>n$; $p_k=\sum_{i=1}^n x_i^k$.
If $n=p$, the statement is true if you require your conditions for all $k$, not just $k\le n$.
Indeed, Newton's identities say that $$ ke_k=\sum_{i=1}^k (-1)^{i-1} p_ie_{k-i} $$ for all $k$ Of course for $k=p$ we shall have an issue with finding $e_p$ from the $p$-th equation. But, if you look at the $p+1$-st equation, you can recover $e_p$ from it unless $p_1=0$, or if $p_1=p_2=\cdots=p_{m-1}=0$ and $p_m\ne0$, then the $m+p$-th equation will give you $e_p$, the $m+p+1$st equation will give you $e_{p+1}$ etc. Thus, you recover all the elementary symmetric polynomials unless $p_k=0$ for all $k$. In this case, Newton's identities tell us that $e_k=0$ for $k$ not divisible by $p$, so the polynomial $$ (x-x_1)(x-x_2)\cdots(x-x_n) $$ is a polynomial in $x^p$. Over a finite field of characteristic $p$, we have $g(x^p)=(h(x))^p$ for some $h(x)$, since every element is a $p$-th power. Thus, there will be repetitions among $x_i$, which is what you do not allow.
Originally, I claimed this for all $n$, but as noted in comments, for $n>p$ one has to be more careful.
In addition, for $n=3=p$, the formula given by @zibadawatimmy here can be made explicit: $(0,1,-1)$ and $(-1+i,i,1+i)$ have the same first three power sums in $\mathbb{F}_9=\mathbb{F}_3[i]/(i^2+1)$. Indeed: \begin{gather} 0+1+(-1)=0=(-1+i)+i+(1+i),\\ 0^2+1^2+(-1)^2=-1=(-1+i)^2+i^2+(1+i)^2,\\ 0^3+1^3+(-1)^3=0=(-1+i)^3+i^3+(1+i)^3. \end{gather} Thus, it is not sufficient to ask for these formulas for $1\le k\le n$.
A counterexample: Let $\mathbf{F}_4 = \{0,1,\alpha,\alpha+1\}$, take $X = \{0,\alpha+1\}$ and $Y = \{1,\alpha\}$. The sums obviously match, and we have
$0^2 + (\alpha+1)^2 = \alpha = 1^2 + \alpha^2$.
Hmmm... If we go to the next higher power sum, though, we have
$0^3 + (\alpha+1)^3 = 1 \ne 0 = 1^3 + \alpha^3$.
Edit: Oops, somehow I missed the requirement to have characteristic more than $2$. To salvage my honor, here's a counterexample in $\mathbb{F}_{p^p}$ with $|X|=|Y|=p$ disjoint:
Take $X = \{0, 1, ..., p-1\}$, so $X$ is the set of solutions to the equation $x^p=x$, and take $Y$ to be the set of solutions to the equation $y^p = y+1$ ($Y$ is a shift of $X$, so the elements of $Y$ are distinct). For $k < p$, the power sums match by Newton's identities, while for $k = p$ we have
$\sum x_i^p = \sum x_i = \sum y_i = \sum (y_i + 1) = \sum y_i^p$.
In fact, by the same argument, the first $k$ such that the power sum differs is $k = 2p-1$.
It's not true. You can easily check it in GAP, for example. Here's a way to get a counterexample in $\mathbb{F}_9$ for $n=4$.
q:=9; n:=4;
FF:=GF(q);
sets:=Tuples(FF,n);;
sets:=Set(List(sets,Set));;
sets:=Filtered(sets,x->Length(x)=n);;
PowerSums := function(vals,m)
return(List([1..m],k->Sum(List(vals,t->t^k))));
end;
S:=First(sets,x->ForAny(sets,y->(not x=y) and PowerSums(x,n)=PowerSums(y,n)));
This returns a four element set, and we can find the other one similarly.
T:=First(sets,y->(not y=S) and PowerSums(S,n)=PowerSums(y,n));
We also get examples for $n=3$ by just changing that one value at the start and rerunning.
Let me prove that if for two sets of $n$ distinct numbers $X=\{x_1,\ldots,x_n\}$ and $Y=\{y_1,\ldots,y_n\}$ the sums of $k$th powers for $k=1,\ldots,2n-1$ are the same, then $X=Y$, irrespectively of the ground field. This shows that the example of @zeb given in his answer (for $n=p$ over a field of characteristic $p$) is, in a sense, optimal.
For that, let us take set $X$ and consider the Newton identities number $n$, \ldots, $2n-1$: \begin{aligned} ne_n&=\sum_{k=1}^n (-1)^{k-1}p_k e_{n-k},\\ 0&=\sum_{k=1}^{n+1} (-1)^{k-1}p_k e_{n+1-k},\\ \ldots \\ 0&=\sum_{k=1}^{2n-1} (-1)^{k-1}p_k e_{2n-1-k},\\ \end{aligned} with the usual convention $e_0=1$, $e_m=0$ for $m>n$. Given the power sums $p_1, p_2, \ldots$, these formulas are $n$ linear equations with $n$ unknowns $e_1, \ldots, e_n$, whose matrix of coefficients is easily seen to have the entry $i,j$ equal to $(-1)^{i+j-1}p_{i+j-2}$, where $p_i$ is the $i$th power sum of $x_1,\ldots,x_n$, including $p_0=n$ (this accounts for the term $ne_n$ in the first equation). Changing the variables $x_i\to -x_i$ and multiplying the matrix by $-1$, we get the matrix whose entry $i,j$ equal to $p_{i+j-2}$. This matrix is, by a direct inspection, equal to $VV^T$, where $V$ is the Vandermonde matrix, so its determinant is $\prod_{i<j}(x_i-x_j)^2$. Thus, for pairwise distinct $x_i$, this system has the only solution for $e_1,\ldots,e_n$. Clearly, the elementary symmetric functions determine the set $X$ uniquely.