If your Prolog has some kind of non backtrackable assignment, like SWI-Prolog 'global' variables, you could implement (beware, simple minded code) in this way:
:- meta_predicate solutions(0, ?).
:- meta_predicate solutions(+, 0, ?).
solutions(G, L) :-
solutions(G, G, L).
solutions(P, G, L) :-
( nb_current(solutions_depth, C) -> true ; C=1 ),
D is C+1,
nb_setval(solutions_depth, D),
atom_concat(solutions_depth_, D, Store),
nb_setval(Store, []),
( G,
nb_getval(Store, T),
nb_setval(Store, [P|T]),
fail
; nb_getval(Store, R)
),
nb_delete(Store),
nb_setval(solutions_depth, C),
reverse(R, L).
Usage of 'global' variables results in more efficient execution WRT the dynamic database (assert/retract), and (in SWI-prolog) can be used even in multithreaded applications.
edit
Thanks to @false comment, moved the cut(s) before reverse/2, and introduced a stack for reentrant calls: for instance
?- solutions(X-Ys,(between(1,3,X),solutions(Y,between(1,5,Y),Ys)),S).
S = [1-[1, 2, 3, 4, 5], 2-[1, 2, 3, 4, 5], 3-[1, 2, 3, 4, 5]].
edit
Here is a variant of solutions/3, building the result list in order, to avoid the final reverse/2 call. Adding results to the list tail is a bit tricky...
solutions(P, G, L) :-
( nb_current(solutions_depth, C) -> true ; C=1 ),
D is C+1,
nb_setval(solutions_depth, D),
atom_concat(solutions_depth_, D, Store),
( G,
( nb_current(Store, U/B) -> B = [P|R], Q = U/R ; Q = [P|T]/T ),
nb_setval(Store, Q),
fail
; ( nb_current(Store, L/[]) -> true ; L = [] )
),
nb_delete(Store),
nb_setval(solutions_depth, C).