As Gianluca said you have much more equations than unknowns, but you may consider to choose a least square fit approach.
Your unknowns are
v = {c1, c2, c3, c4, c5, c6, c7, c8, c9};
and your equations
eqs = List[{uj == aj*c1 + bj*c2 + cj*c3 + ui,
vj == aj*c4 + bj*c5 + cj*c6 + vi,
tj == aj*c7 + bj*c8 + cj*c9 + ti,
uk == ak*c1 + bk*c2 + ck*c3 + ui,
vk == ak*c4 + bk*c5 + ck*c6 + vi,
tk == ak*c7 + bk*c8 + ck*c9 + ti,
ul == al*c1 + bl*c2 + cl*c3 + ui,
vl == al*c4 + bl*c5 + cl*c6 + vi,
tl == al*c7 + bl*c8 + cl*c9 + ti,
um == am*c1 + bm*c2 + cm*c3 + ui,
vm == am*c4 + bm*c5 + cm*c6 + vi,
tm == am*c7 + bm*c8 + cm*c9 + ti,
un == an*c1 + bn*c2 + cn*c3 + ui,
vn == an*c4 + bn*c5 + cn*c6 + vi,
tn == an*c7 + bn*c8 + cn*c9 + ti,
uo == ao*c1 + bo*c2 + co*c3 + ui,
vo == ao*c4 + bo*c5 + co*c6 + vi,
to == ao*c7 + bo*c8 + co*c9 + ti,
up == ap*c1 + bp*c2 + cp*c3 + ui,
vp == ap*c4 + bp*c5 + cp*c6 + vi,
tp == ap*c7 + bp*c8 + cp*c9 + ti}][[1]];
Now make a function to transform your equations by bringing the additive constants to the left hand sides
trans[e_] := # - Select[e[[2]], FreeQ[#, Times] &] & /@ e
and get the modified equations ( you could have written your equations directly in this form )
eqs1 = trans /@ eqs
Make a function to extract the coefficients of your unknowns in matrix-form
tomat[e_] := Coefficient[e[[2]], #] & /@ v
and get that matrix
mat = tomat /@ eqs1;
% // MatrixForm
the known components of yor "result"-vector are
vec = #[[1]] & /@ eqs1
This vector and the matrix indeed yield your former system of equations
Thread[vec == mat.v]
For the least square fit define a deviation vector dev
dev = mat . v - vec
and find v so that dev.dev (as fuction of v = {c1, c2, ... , c9 } ) attains a minimum.
If I did not make a mistake in the caclulations this solution v should be given by
v = Inverse[ Transpose[ mat ] . mat ]. Transpose[ mat ] . vec
You can do this symbolically, but that yields outputs of enormous lenghts, It seems to be better to insert numbers for your coeffcients.
