R, 69 bytes

function(n,b=n+1){while(sum(T)-n)T=((F=F+1)%/%(b^(0:3))%%b)^2;T[T>0]}

Try it online!

Counts up from 1, converting each number to base-n digits, least-significant first, and returning the first set for which the sum of squares of the digits equals n.

(previous version, with 4 bytes saved thanks to Robin Ryder):
# R, 75 71 bytes

function(n,x=expand.grid(a<-0:n,a,a,a)^2,y=x[rowSums(x)==n,][1,])y[y>0]

Try it online!

Uses a variation of Robin Ryder's approach, here constructing all combinations of four squares in one step using the built-in R function expand.grid.

R, 67 bytes

function(n,b=n+1){while(sum(T)-n)T=((F=F+1)%/%b^(0:3)%%b)^2;T[T>0]}

Try it online!

Counts up from 1, converting each number to base-n digits, least-significant first, and returning the first set for which the sum of squares of the digits equals n.

(previous version, with 4 bytes saved thanks to Robin Ryder):
# R, 75 71 bytes

function(n,x=expand.grid(a<-0:n,a,a,a)^2,y=x[rowSums(x)==n,][1,])y[y>0]

Try it online!

Uses a variation of Robin Ryder's approach, here constructing all combinations of four squares in one step using the built-in R function expand.grid.

R, 75 bytes

function(n,x=expand.grid(rep(list(0:n),4))^2,y=x[rowSums(x)==n,][1,])y[y>0]

Try it online!

Uses a variation of Robin Ryder's approach, here constructing all combinations of four squares in one step using the built-in R function expand.grid.

R, 69 bytes

function(n,b=n+1){while(sum(T)-n)T=((F=F+1)%/%(b^(0:3))%%b)^2;T[T>0]}

Try it online!

Counts up from 1, converting each number to base-n digits, least-significant first, and returning the first set for which the sum of squares of the digits equals n.

(previous version, with 4 bytes saved thanks to Robin Ryder):
# R, 75 71 bytes

function(n,x=expand.grid(a<-0:n,a,a,a)^2,y=x[rowSums(x)==n,][1,])y[y>0]

Try it online!

Uses a variation of Robin Ryder's approach, here constructing all combinations of four squares in one step using the built-in R function expand.grid.