Using the solution of a linear system, and splitting a matrix

edit

Hi there, I'm fairly new to Sage and Python, so I'm getting into basic problems here that I'd be happy if you could help me out.
Here it is, actually, here they are: I'm generating a set of equalities and solving them with solve. A simple example:

import numpy as n;
m=4;
s = list(var('s_%d' % int(i)) for i in range(m));
eqns=[s_0+s_1==1,s_2-s_3==0];
sol=solve(eqns,s,solution_dict=True)[0]

This gives the solutions:

{s_1: -r2 + 1, s_0: r2, s_3: r1, s_2: r1}

My first question is, how do I create a matrix with the solutions? Say, something like:

M=m.zeros((2,2));
for i in range(2):
for j in range(2):
M[int(i),int(j)] = sol[s[i]]+ sol[s[j+2]]

This is giving me the error: "TypeError: unable to simplify to float approximation"

My second question, would be, given the array M, how do I split it as r1 times a matrix, plus r2 times another matrix, plus a constant matrix? In the above example,

M= [[1,r1+r2],[r1-r2+1,r1-r2+1]]= r1 [[0,1],[1,1]] + r2 [[0,1],[-1,-1]]+ [[1,0],[0,0]]

I'm interested in the matrices multiplying the still unknown coefficients. Maybe I should add that the number of equations in the problem I'm solving is much bigger than in this simple example, and therefore I cannot find this matrix decomposition by simply looking at it.

Thanks for the help!