In three dimensions, a $3\times 3$ rotation matrix can perform any conceivable rotation you might want upon a vector.
But a vector by itself has no "location" about which to rotate. It merely changes direction. Even a large collection of vectors still has no "location."
The fact that you want to rotate a model that has a position means you are dealing with spatial coordinates, not merely vectors.
One way to deal with spatial coordinates, of course, is to nominate one point in space as the "origin" and identify every point $P$ with the vector from the origin to $P.$
When you do that, rotations still change the directions of the vectors;
since we still identify each point by its vector from the origin,
rotating that vector rotates that point around an axis through the origin.
This means the origin, which is represented by a zero vector (all coordinates zero), never moves when you apply a $3\times 3$ rotation matrix in three dimensions.
In fact, no matter what $3\times 3$ matrix you multiply with the zero vector you still get back the zero vector, even if the matrix does not represent a true rotation.
If you want to rotate a model around its own center which is not the origin,
then unless you happen to pick the one rotation axis that also goes through the origin, your rotation will move the origin.
Hence it is impossible to do with a $3\times 3$ matrix.
What you can do, if you want, is to construct a rotation $R$ by the desired angle around an axis in the desired orientation (parallel to the axis you really want).
Then if the point $C$ is the center of your object, compute $R(C).$
Now find the vector $v = C - R(C).$
Then you can perform the desired rotation of any point $P$ around your object's own center by applying the rotation $R$ and the the translation $T_v$ (adding the vector $v$ to the point's coordinates):
$$ T_v(R(P)). $$
Now it's just one translation in addition to the rotation.
The "axis-angle rotation" matrix you showed, which is a $4\times 4$ matrix,
is actually more than just a rotation matrix, however.
The way it works is that any point $P$ is represented by a four-coordinate vector,
x_P\\ y_P\\ z_P\\ 1,
which you can rotate by multiplying by a matrix in the form you gave;
but in order to follow the rotation by a translation, you just need to replace the zeros in the right-hand column of your matrix by the coordinates of the translation vector.
You can also just take the upper left $3\times 3$ submatrix of your "axis-angle" matrix and multiply it by the three-coordinate vector from the origin to each point,
but then it's just a rotation around the origin like any other
$3\times 3$ matrix rotation in three dimensions.
Whether this avoids "gimbal lock" depends on what you want to do with it.
If you're actually using it to track the rotation of a three-axis gimbal,
you can still "lock" the gimbal by putting it through some sequence of axis-angle rotations.
If you have a free-floating object in space, not an actual gimbal,
then the way in which "gimbal lock" might apply is if you used a set of Euler angles to represent every orientation of the object relative to its original orientation.
Some sequence of orientations will force you to suddenly "flip" one of your Euler angles, which is the kind of manipulation that "locks" a real-life gimbal.
That is, the Euler angles are a kind of fictional gimbal.
But as long as you don't have a real gimbal or a fictional gimbal,
you can't "lock" it.
The difference between matrices and quaternions is an entire question of its own.
But when they are applied to three-dimensional rotations, each in their own usual way, each is just a different way to represent a rotation.
And each has its own way of calculating the results of the rotation.
There is plenty written about either of these techniques that you can study.