Prolog Sudoku solver taking too long
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.