coloring a $4\times 4$ square with $4$ colors
As requested, using the polynomial $$f(t):= t ( t-1 ) ( {t}^{14}-23\,{t}^{13}+253\,{t}^{12}-1762 \,{t}^{11}+8675\,{t}^{10}-31939\,{t}^{9}+90723\,{t}^{8}-202160\,{t}^{7 }+355622\,{t}^{6}-492434\,{t}^{5}+529770\,{t}^{4}-430857\,{t}^{3}+ 251492\,{t}^{2}-94782\,t+17493 )$$ and calculating $f(4)-4f(3)+6f(2)$, the answer is 5969496 by Jonas' observation. The 4 is the four different sets of three of the colours and the 6 is the pairs of two colours.
The calculations for calculating the chromatic polynomial of the $4\times 4$ grid are too involved to post here, though. Hopefully I can show the ideas behind a $3\times 3$ and you can imagine the work required to do the work for the question by hand.
My chromatic polynomial reduction rules are given on pages 24 and 25 of my course notes
Using deletion-contraction on a central edge (1:2- 2:2) you get that $g(t) = g_1(t) - g_2(t)$ where $G_1$ and $G_2$ are the following graphs:
We can then use deletion-contraction again on the opposite edge in both $G_1$ and $G_2$ and get $f(t) = h_1(t) - 2 h_2(t) + h_3(t)$ where the graphs are as follows:
From here we can use complete intersection and maximum valency for $H_2$ and $H_3$, respectively so that $h_2(t) = (t-2)^2 \times p(C_6,t) = t \left( t-1 \right) \left( {t}^{4}-5\,{t}^{3}+10\,{t}^{2}-10\,t+5 \right) \left( t-2 \right) ^{2} $ and $h_3(t) = t \times p(P_3 \cup P_3,t-1) = t \left( t-1 \right) ^{2} \left( t-2 \right) ^{4}$, using the standard results for polynomials of cycles, paths and unions.
One more deletion-contraction gives $h_1 = t \left( t-1 \right) \left( {t}^{7}-9\,{t}^{6}+36\,{t}^{5}-84\,{t}^{ 4}+126\,{t}^{3}-124\,{t}^{2}+76\,t-23 \right) $ and putting them together gives us $$g(t) = t \left( t-1 \right) \left( {t}^{7}-11\,{t}^{6}+55\,{t}^{5}-161\,{t} ^{4}+298\,{t}^{3}-350\,{t}^{2}+244\,t-79 \right) $$
Note that $g(2)=2$ since the graph is bipartite, $g(3) =246$ and $g(4) =9612$, so that there are 8640 colourings using exactly 4 colours.