A: matrix(
 [1 -l -s*l_s,1-l -s*l_s - Nsl_S,M*(l_s + Nl_S)], 
 [-s*l_S,-Nsl_S,M*Nl_S], 
 [-s*l_D,-NSl_D,1-delta + M*Nl_D + g_D]
);

invert(A);

b: matrix(
 [b1], 
 [b2], 
 [b3]
);

x : invert(A).b;

x[1];

opt_cond : f/beta = 2*v_St + m*v_Dt + v_st;

subst(v_st=x[1],opt_cond);

determinant(A);

adjoint(A).b;

test : matrix([1,2],[3,4]);