solving system of ODE and exporting results into matrix form

> 1) format results in form "diff(t1(tau),tau)
> =A*t1(tau)+B*t2(tau)+C*q1(tau)"
> 
> 2) substitute variables in results "lambda/(c*ro*delta^2)" with "b" and
> "1/(c*ro*delta)" with "d"
> 
> 3) represent system of ODE in "dt_dtau=F*t+G*u" matrix form, where:
> dt_dtau:transpose(matrix([diff(t1(tau),tau),diff(t2(tau),tau),diff(t3(tau),tau)]));
> t:transpose(matrix([t1(tau),t2(tau),t3(tau)]));
> u:transpose(matrix([q1(tau),q2(tau)]));
> 
> Thanks in advance.
> 

coefmatrix() does your job. "Automation" of introducing the 
abbreviations c, b for arbitrary equations should be difficult.

wxMaxima 0.8.1 http://wxmaxima.sourceforge.net
Maxima 5.17.1 http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (aka GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
/*DEs*/
(%i1) 
eq1:q1(tau)=c*ro*S*delta/2*diff(t1(tau),tau)+lambda*S*(t1(tau)-t2(tau))/delta;
eq2:lambda*S*(t1(tau)-t2(tau))/delta=c*ro*S*delta*diff(t2(tau),tau)+lambda*S*(t2(tau)-t3(tau))/delta;
eq3:lambda*S*(t2(tau)-t3(tau))/delta=c*ro*S*delta/2*diff(t3(tau),tau)+q2(tau);
(%o1) 
q1(tau)=((t1(tau)-t2(tau))*S*lambda)/delta+(c*delta*ro*('diff(t1(tau),tau,1))*S)/2
(%o2) 
((t1(tau)-t2(tau))*S*lambda)/delta=((t2(tau)-t3(tau))*S*lambda)/delta+c*delta*ro*('diff(t2(tau),tau,1))*S
(%o3) 
((t2(tau)-t3(tau))*S*lambda)/delta=(c*delta*ro*('diff(t3(tau),tau,1))*S)/2+q2(tau)
/*sols*/
(%i4) s1:solve(eq1, diff(t1(tau),tau));
s2:solve(eq2, diff(t2(tau),tau));
s3:solve(eq3, diff(t3(tau),tau));
(%o4) 
['diff(t1(tau),tau,1)=((2*t2(tau)-2*t1(tau))*S*lambda+2*delta*q1(tau))/(c*delta^2*ro*S)]
(%o5) 
['diff(t2(tau),tau,1)=((t3(tau)-2*t2(tau)+t1(tau))*lambda)/(c*delta^2*ro)]
(%o6) 
['diff(t3(tau),tau,1)=-((2*t3(tau)-2*t2(tau))*S*lambda+2*delta*q2(tau))/(c*delta^2*ro*S)]
/*vecs*/
(%i7) dtdtau: map(lambda([e],part(e,1,1)),[s1,s2,s3]);
t: map(lambda([e], part(e,1)), dtdtau);
u: [q1(tau),q2(tau)];
(%o7) ['diff(t1(tau),tau,1),'diff(t2(tau),tau,1),'diff(t3(tau),tau,1)]
(%o8) [t1(tau),t2(tau),t3(tau)]
(%o9) [q1(tau),q2(tau)]
(%i10) eqset: map(
     lambda([e], lambda([L], L[2]-L[1])(args(part(e,1)))),
     [s1,s2,s3]
);
(%o10) 
[((2*t2(tau)-2*t1(tau))*S*lambda+2*delta*q1(tau))/(c*delta^2*ro*S)-'diff(t1(tau),tau,1),((t3(tau)-2*t2(tau)+t1(tau))*lambda)/(c*delta^2*ro)-'diff(t2(tau),tau,1),-((2*t3(tau)-2*t2(tau))*S*lambda+2*delta*q2(tau))/(c*delta^2*ro*S)-'diff(t3(tau),tau,1)]
/*mats*/
(%i11) Fmat: coefmatrix(eqset,args(t));
(%o11) 
matrix([-(2*lambda)/(c*delta^2*ro),(2*lambda)/(c*delta^2*ro),0],[lambda/(c*delta^2*ro),-(2*lambda)/(c*delta^2*ro),lambda/(c*delta^2*ro)],[0,(2*lambda)/(c*delta^2*ro),-(2*lambda)/(c*delta^2*ro)])
(%i12) Gmat: coefmatrix(eqset, u);
(%o12) matrix([2/(c*delta*ro*S),0],[0,0],[0,-2/(c*delta*ro*S)])
(%i13) dtdtaumat: transpose(dtdtau);
tmat: transpose(t);
umat: transpose(u);
(%o13) 
matrix(['diff(t1(tau),tau,1)],['diff(t2(tau),tau,1)],['diff(t3(tau),tau,1)])
(%o14) matrix([t1(tau)],[t2(tau)],[t3(tau)])
(%o15) matrix([q1(tau)],[q2(tau)])
(%i16) Fmat . tmat + Gmat . umat - dtdtaumat;
(%o16) 
matrix([(2*t2(tau)*lambda)/(c*delta^2*ro)-(2*t1(tau)*lambda)/(c*delta^2*ro)+(2*q1(tau))/(c*delta*ro*S)-'diff(t1(tau),tau,1)],[(t3(tau)*lambda)/(c*delta^2*ro)-(2*t2(tau)*lambda)/(c*delta^2*ro)+(t1(tau)*lambda)/(c*delta^2*ro)-'diff(t2(tau),tau,1)],[-(2*t3(tau)*lambda)/(c*delta^2*ro)+(2*t2(tau)*lambda)/(c*delta^2*ro)-(2*q2(tau))/(c*delta*ro*S)-'diff(t3(tau),tau,1)])
/*did this work?*/
(%i17) % - transpose(eqset);
ratsimp(%);
(%o17) 
matrix([-((2*t2(tau)-2*t1(tau))*S*lambda+2*delta*q1(tau))/(c*delta^2*ro*S)+(2*t2(tau)*lambda)/(c*delta^2*ro)-(2*t1(tau)*lambda)/(c*delta^2*ro)+(2*q1(tau))/(c*delta*ro*S)],[-((t3(tau)-2*t2(tau)+t1(tau))*lambda)/(c*delta^2*ro)+(t3(tau)*lambda)/(c*delta^2*ro)-(2*t2(tau)*lambda)/(c*delta^2*ro)+(t1(tau)*lambda)/(c*delta^2*ro)],[((2*t3(tau)-2*t2(tau))*S*lambda+2*delta*q2(tau))/(c*delta^2*ro*S)-(2*t3(tau)*lambda)/(c*delta^2*ro)+(2*t2(tau)*lambda)/(c*delta^2*ro)-(2*q2(tau))/(c*delta*ro*S)])
(%o18) matrix([0],[0],[0])
/*abbrevs*/
(%i19) Fmat[1,1]=-c;
replacelambda: first(solve(%, lambda));
Gmat[1,1]=2*b;
replacero: first(solve(%, ro));
(%o19) -(2*lambda)/(c*delta^2*ro)=-c
(%o20) lambda=(c^2*delta^2*ro)/2
(%o21) 2/(c*delta*ro*S)=2*b
(%o22) ro=1/(b*c*delta*S)
(%i23) F: subst(replacelambda, Fmat);
(%o23) matrix([-c,c,0],[c/2,-c,c/2],[0,c,-c])
(%i24) G: subst(replacero, Gmat);
(%o24) matrix([2*b,0],[0,0],[0,-2*b])
/*Equation*/
(%i25) box(dtdtaumat = dot(F, tmat) + dot(G, umat));
(%o25) 
matrix(['diff(t1(tau),tau,1)],['diff(t2(tau),tau,1)],['diff(t3(tau),tau,1)])=dot(matrix([-c,c,0],[c/2,-c,c/2],[0,c,-c]),matrix([t1(tau)],[t2(tau)],[t3(tau)]))+dot(matrix([2*b,0],[0,0],[0,-2*b]),matrix([q1(tau)],[q2(tau)]))