restart;with(LinearAlgebra):
ep[0]:=0;
for k from 1 to 3 do
   qp[k]:=((q[k]-ep[k-1])-v)+e[k];
   ep[k]:=e[k]*(q[k+1]/qp[k]);
end do;
qp[4]:=(q[4]-ep[3])-v;
As:=Matrix(4,4,[[qp[1],-qp[1],0,0],[-ep[1],qp[2]+ep[1],-qp[2],0],
                [0,-ep[2],qp[3]+ep[2],-qp[3]],[0,0,-ep[3],qp[4]+ep[3]]]);
D1:=DiagonalMatrix([1,ep[1]/e[1], ep[1]/e[1]*ep[2]/e[2],
                    ep[1]/e[1]*ep[2]/e[2]*ep[3]/e[3]]);
C:=MatrixInverse(D1).As.D1;
X:=Matrix(4,4,[[q[1],0,0,0],[-e[1],q[2],0,0],[0,-e[2],q[3],0],
               [0,0,-e[3],q[4]]]);
Y:=Matrix(4,4,[[1,-1,0,0],[0,1,-1,0],
               [0,0,1,-1],[0,0,0,1]]);
Bv:=Y.X-v-C;
simplify(%);