Well, copying my answer from SO.
Think about the number of outputs of subset/4, if every cell of Solution is unbound in the beginning. There are 44*4 = 4,294,967,296 outputs. Now you might know why your code is so slow.
You should call valid/1 inside of subset/2. Test if the puzzle is still valid after inserting one number:
subset(_,[],_).
subset(RCS,[H|T],List):-
member(H,List),
valid(RCS),
subset(RCS,T,List).
% Like \+member(H,T) but unbound cells are ignored
nonmember(_, []).
nonmember(H1, [H2|T]) :-
(var(H2); H1 \= H2),
nonmember(H1, T).
different([]).
different([H|T]):-
(var(H); nonmember(H,T)),
different(T),!.
valid([]).
valid([Head|Tail]) :-
different(Head),
valid(Tail),!.
sudoku(Puzzle) :-
Puzzle = [S11, S12, S13, S14,
S21, S22, S23, S24,
S31, S32, S33, S34,
S41, S42, S43, S44],
Row1 = [S11, S12, S13, S14],
Row2 = [S21, S22, S23, S24],
Row3 = [S31, S32, S33, S34],
Row4 = [S41, S42, S43, S44],
Col1 = [S11, S21, S31, S41],
Col2 = [S12, S22, S32, S42],
Col3 = [S13, S23, S33, S43],
Col4 = [S14, S24, S34, S44],
Square1 = [S11, S12, S21, S22],
Square2 = [S13, S14, S23, S24],
Square3 = [S31, S32, S41, S42],
Square4 = [S33, S34, S43, S44],
RCS = [Row1, Row2, Row3, Row4,
Col1, Col2, Col3, Col4,
Square1, Square2, Square3, Square4],
subset(RCS, Puzzle, [1,2,3,4]).
Next optimization would be not to test all rows, columns and squares after every insert, but only the affected ones (3 instead of 12):
subset([], _).
subset([(P,RCS) | T], List):-
member(P, List),
valid(RCS),
subset(T, List).
% Like \+member(H,T) but unbound cells are ignored
nonmember(_, []).
nonmember(H1, [H2|T]) :-
(var(H2); H1 \= H2),
nonmember(H1, T).
different([]).
different([H|T]):-
(var(H); nonmember(H,T)),
different(T),!.
valid([]).
valid([Head|Tail]) :-
different(Head),
valid(Tail),!.
sudoku(Puzzle) :-
Puzzle = [S11, S12, S13, S14,
S21, S22, S23, S24,
S31, S32, S33, S34,
S41, S42, S43, S44],
Row1 = [S11, S12, S13, S14],
Row2 = [S21, S22, S23, S24],
Row3 = [S31, S32, S33, S34],
Row4 = [S41, S42, S43, S44],
Col1 = [S11, S21, S31, S41],
Col2 = [S12, S22, S32, S42],
Col3 = [S13, S23, S33, S43],
Col4 = [S14, S24, S34, S44],
Square1 = [S11, S12, S21, S22],
Square2 = [S13, S14, S23, S24],
Square3 = [S31, S32, S41, S42],
Square4 = [S33, S34, S43, S44],
RCS11 = (S11, [Row1, Col1, Square1]),
RCS12 = (S12, [Row1, Col2, Square1]),
RCS13 = (S13, [Row1, Col3, Square2]),
RCS14 = (S14, [Row1, Col4, Square2]),
RCS21 = (S21, [Row2, Col1, Square1]),
RCS22 = (S22, [Row2, Col2, Square1]),
RCS23 = (S23, [Row2, Col3, Square2]),
RCS24 = (S24, [Row2, Col4, Square2]),
RCS31 = (S31, [Row3, Col1, Square3]),
RCS32 = (S32, [Row3, Col2, Square3]),
RCS33 = (S33, [Row3, Col3, Square4]),
RCS34 = (S34, [Row3, Col4, Square4]),
RCS41 = (S41, [Row4, Col1, Square3]),
RCS42 = (S42, [Row4, Col2, Square3]),
RCS43 = (S43, [Row4, Col3, Square4]),
RCS44 = (S44, [Row4, Col4, Square4]),
RCS = [RCS11, RCS12, RCS13, RCS14,
RCS21, RCS22, RCS23, RCS24,
RCS31, RCS32, RCS33, RCS34,
RCS41, RCS42, RCS43, RCS44],
subset(RCS, [1,2,3,4]).
Now you should be able to increase the puzzle size to 9×9 and still have a reasonably fast solver.