A = np.matrix([[deltaItot/(2*A1), -Lf/tf*1/(2*A1), -1],
               [-Lf/tf*1/(2*A2), deltaIItot/(2*A2), -1],
               [2*A1, 2*A2, 0]])
b = np.matrix([ [ -np.sum(np.multiply(qbI,deltaI))/(2*A1)],
              [ - np.sum(np.multiply(qbII,deltaII))/(2*A2)],
               [ - (np.sum(np.multiply(qbI[:20],LI[:20])*R) \
         + np.sum(np.multiply(qbI[25:],LI[25:])*R) + \
         np.sum(np.multiply(qbI[20:25],LI[20:25])*(R-hf)) \
         + np.sum(np.multiply(qbII[:16],LII[:16])*R) + \
         np.sum(np.multiply(qbII[16:],LII[16:])*(R-hf)))]])
c1 = np.linalg.solve(A,b)
print c1