I derived P2 before. P3 should be similar.
You already know the dimension, which is good. We also need to know the position of nodes. For P2, we have 4 at vertices, and 6 at edges.
Then the rest is to determine the coefficients of an interpolating polynomial.
Here is the general idea.
The 10 monomials $m_i$ are 1, x, y, z, x*y, ... etc. For each term we want to determine its coefficient $c_i$,
So that interpolating polynomial $q = \sum_i m_i c_i$ evaluates exact node position given x, y, and z.
Solving the system will give us $c_i$, and substitute them back to $q$, regrouping the terms will produce the basis in terms of x,y,and z.
Here is some maple code to help with the derivation:
# the 10 monomials
# We are going to determine ci
q:=`+`(seq(c||i * m[i],i=1..10));
# generate node positions
# Evaluate at node coordinates, and require them to be equal to node positions qi.
# Solve for coefficients ci, and substitute back to q
# Find the coefficient of each qi (the node position)
Then the P2 basis from above code:
[2*x^2+4*x*y+4*x*z+2*y^2+4*y*z+2*z^2-3*x-3*y-3*z+1, 2*x^2-x, 2*y^2-y, 2*z^2-z, -4*x^2-4*x*y-4*x*z+4*x, -4*x
*y-4*y^2-4*y*z+4*y, -4*x*z-4*y*z-4*z^2+4*z, 4*x*y, 4*z*x, 4*y*z]