Equation not fully solved for

I think I've worked out the details. I had to do some stuff by hand, and Maxima was mostly just checking that I did it correctly. If the goal is just to get to a solution, I guess that's okay.

Here's my script. You can execute this via: maxima --batch=foo.mac where foo.mac is the name of the script.

/* adapted from: https://stackoverflow.com/questions/69700162/equation-not-fully-solved
 */

load("vect");

/* I don't want to reorder arguments of cross product. */
remrule ("~", ?\~rule4);

/* I do want to flatten noncommutative multiplication.
 * (It appears that's disabled by vect; ordinarily it happens by default.)
 */
dotassoc: true $

declare(P, nonscalar);
declare([v1,w1,V1,W1,r1], nonscalar);
declare([v2,w2,V2,W2,r2], nonscalar);
declare([m1,I1,m2,I2], scalar);

/* Revolute Constraint */
constraint: v2 - (r2~w2) - (v1 - (r1~w1)) = 0$

/* Velocities after impulse P */
eq1: v1 = V1 - P/m1$
eq2: w1 = W1 - (r1~P) / I1$
eq3: v2 = V2 + P/m2$
eq4: w2 = W2 + (r2~P) / I2$

eq: subst([eq1,eq2,eq3,eq4], constraint);

A(a) := matrix([0, -a[3], a[2]], [a[3], 0, -a[1]], [-a[2], a[1], 0]);

eqa: ev (eq, (u~v) := 'A(u).v);

matchdeclare (pp, lambda ([e], not freeof(P, e)));
matchdeclare (qq, lambda ([e], freeof(P, e)));
defrule (rp, pp + qq = 0, pp = -qq);

eqa1: expand (eqa);
eqa2: apply1 (eqa1, rp);

matchdeclare (aa, lambda ([e], matrixp(e) or listp(e)));
tellsimpafter (I(aa), ident (length (aa)));

matchdeclare ([aa, bb], all);
tellsimpafter (I(aa) . bb, bb);
tellsimpafter (aa . I(bb), aa);

M: -(1/I2)*'A(r2).'A(r2) - (1/I1)*'A(r1).'A(r1) + (1/m2)*I(P) + (1/m1)*I(P);
N: 'A(r2).W2 - 'A(r1).W1 - V2 + V1;

eqa2_factored: M . P = N;

expand (eqa2_factored);
?resimplify (%);
if % # eqa2 then error ("tried to factor eqa2, but something went wrong.");

solution: P = M^^-1 . N;

/* EXAMPLE: */

I1: 20 $
I2: 3 $
m1: 100 $
m2: 12 $
V1: [17, 19, -23] $
V2: [-5, -3, 11] $
W1: [8, 4, 14] $
W2: [-6, -16, 24 ] $
r1: [1/2, 2/3, 3/4] $
r2: [5, 7, 3] $

/* note various subterfuges to ensure evaluation with stated values */

example_M: ev (subst (I(P) = ident(3), M), nouns);
example_N: ev (N, nouns);

example_P: example_M^^-1 . example_N;

subst (P = example_P, ev (eqa2, eval, nouns));

if lhs(%) = rhs(%)
    then print ("TEST PASSED: lhs = rhs")
    else error ("TEST FAILED: lhs # rhs");

If you need to evaluate P for different parameters r1, r2, etc., my advice is to evaluate matrix M and vector N with whatever values you want to plug in, and then solve the equation P = M^^-1 . N. An explicit solution is probably going to be pretty messy.

Following Robert Dodier's advice, I broke up all the vectors and solved for P[1] and P[2] individually. I've got something that gives me an answer but now how can I get it into nice vector form? Here it is:

load("vect");
declare(".", commutative);

declare(P, nonscalar);
declare([v1,V1,r1], nonscalar);
declare([v2,V2,r2], nonscalar);
declare([w1,W1,m1,I1], scalar);
declare([w2,W2,m2,I2], scalar);

cross_scalar_vector(s,v) := [-s*v[2], s*v[1]]$
/* Revolute Constraint on Linear Velocity */
constraint: v2 + cross_scalar_vector(w2,r2) - (v1 + cross_scalar_vector(w1,r1)) = [0,0]$

/* Sub in velocities after impulse P. */
post_velocities: [
    v1 = V1 - P/m1,
    w1 = W1 - (r1~P) / I1,
    v2 = V2 + P/m2,
    w2 = W2 + (r2~P) / I2
]$
constraint: subst(post_velocities, constraint)$

/* Break up the remaining vectors for solve. */
vectors: [
    P = [P[1], P[2]],
    V1 = [V1[1], V1[2]],
    r1 = [r1[1], r1[2]],
    V2 = [V2[1], V2[2]],
    r2 = [r2[1], r2[2]]
]$
constraint: subst(vectors, constraint)$

/* Break up vector constraint into x and y constraint for solve. */
xconstraint: lhs(constraint)[1] = 0$
yconstraint: lhs(constraint)[2] = 0$

/* Not sure why we need to do this again? */
xconstraint: subst(vectors, xconstraint)$
yconstraint: subst(vectors, yconstraint)$

/* Expand cross products for solve. */
xconstraint: express(xconstraint)$
yconstraint: express(yconstraint)$

solve([xconstraint,yconstraint], [P[1],P[2]]);