Let $(X_1,X_2,\ldots,X_n)$ be a random sample of size $n$ drawn from a population having pmf $P(X=j)=\frac{1}{N}\mathbf1_{j\in \{1,2,\ldots,N\}}$. What is the joint pmf of $X_{(1)}$ and $X_{(n)}$ ?
Let $U=X_{(1)}$ and $V=X_{(n)}$. Then for $u<v$, the joint distribution function of $(U,V)$ is
$F_{U,V}(u,v)=P(U\leqslant u,V\leqslant v)$
$\qquad\qquad\quad=P(V\leqslant v)-P(U>u,V\leqslant v)$
$\qquad\qquad\quad=P(X_1,X_2,\ldots,X_n\leqslant v)-P(u<X_1,X_2,\ldots,X_n\leqslant v)$
$\qquad\qquad\quad=(F(v))^n-(F(v)-F(u))^n$, where $F$ is the population CDF.
Then, $P(U=u,V=v)=F_{U,V}(u,v)-F_{U,V}(u-1,v)-F_{U,V}(u,v-1)+F_{U,V}(u-1,v-1)$ where $u<v$. I get $F(x)=\frac{x}{N}$ for $x=1,2,\ldots,N$. So I finally end up with
$$P(U=u,V=v)=\frac{1}{N^n}((v-u+1)^n+(v-u-1)^n-2(v-u)^n)\mathbf1_{{u<v}{\{1,2,\ldots,N\}}}$$
Is this a proper way to proceed? I don't seem to get the correct marginal distributions from this.
Your work is almost completely correct. The only issue is that you computed $P(U=u,V=v)$ for all $u<v$, but left out the case $u=v$. This needs to be handled separately.
This does give the correct marginal distributions. When you sum over $u$, you get two telescoping series, which results in $P(V=v)=(v^n-(v-1)^n)/N^n$.