In the ring $R =\mathbb Z_{p^k}$ we can give to each element $g \in R$ two "coordinates" $a$ and $b$ as follows: $b = g \;\operatorname{mod} \; p$ and $a = (g-b)/p$. So we have in a unique fashion that $g = ap+b$, with $0 \leq a \lt p^{k-1}$ and $0 \leq b \lt p$ Consequently one can define a function crd in GAP :
crd := function ( g )
return [ (Int( g ) - Int( g ) mod p) / p, Int( g ) mod p ];
end;
Let $\mathbb Z^*_{p^k}$ the group of units of $R$, characterized by the property that their $b$-coordinate $=1$. On this group we define a homomorpism by $h:G \rightarrow G: g \mapsto g^{p-1}$. The set of equations we want to solve is thus to be found in the kernel of this homomorphism. As an example we take $p=7$ and $k=5$ giving the following result in GAP:
gap> p := 7;; k := 5;; n := p^k;; R := ZmodnZ(n);; G := Units(R);;
gap> hom := GroupHomomorphismByFunction(G,G,g->g^(p-1));;
gap> K := Kernel(hom);;
gap> List(K, crd);
[ [ 0, 1 ], [ 2207, 4 ], [ 193, 2 ], [ 2207, 5 ], [ 2400, 6 ], [ 193, 3 ] ]
From this and every other example I verified, I obtained the same result: for every $0 \lt b \lt p$ there is exactly one $0 \leq a \lt p^{k-1}$ such that $(ap + b)^{p-1} = 1 \mod p^k$, but I can't figure out how to determine this $a$.
The easy way is to use Hensel's lemma to lift the result.
The polynomial $f(x)=x^{p-1}-1$ has simple roots at each nonzero $x\in\mathbb{Z}/p$, since $f'(x)=(p-1)x^{p-2}\neq 0$. So starting at $b_0=b$ we compute the lift by a single Newton-Raphson iteration $$ b_{k+1}\equiv b_k-\frac{f(b_k)}{f'(b_k)}\pmod{p^{k+1}}. $$ and your $a$ is just $(b_k-b)/p$, the representative $b_k\in\{1,2,\dots,p^k-1\}$.