Use linear algebra methods of SciPy to solve the three simultaneous equations

My question is that how to write theese equations in array and solve?

from scipy import linalg
import numpy as np
import matplotlib.pyplot as plt

x = np.array[-23,1100,2300],[2300,1500,550],[550,1600,]

I tried to write in the array above, but I couldn't figure out how to replace 'In' and 'Vs2' in the question. Can you help me solve the question?

You want to solve these equations for several voltages, which suggests the use of a for-loop. For clarity, it's usually better to use identifiers for values, thus for instance, R1 rather than 1100. Put the R1 in formulae and let the computer do the simple arithmetic for you.

You may be thinking of using the linalg solve function since you need to solve a square matrix of order three. The unknowns are the currents. Therefore, do the algebra so that you have expressions for the coefficients of the matrix, and for the right side of the equation, in terms of resistances and voltages.

For the matrix (as indicated in the documentation at https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html#scipy.linalg.solve),

a = np.array([[f1(Rs, Vs), f2(Rs, Vs), f3(Rs, Vs)], [...], [...]])

For the vector on the right side,

b = np.array([f4(Rs, Vs), f5(Rs,Vs), f6(Rs, Vs)])

Then currents = solve(a, b)

Notice that f1, f2, etc are those functions that you have to calculate algebraically.

Now put this code in a loop, more or less like this:

for vs2 in [10,15,20,25]:
    currents = solve(a, b)

Because you've got the resistances and vs2's in your algebraic expressions you'll get the corresponding currents. You'll need to collect the currents corresponding to voltages for plotting.

Addition: Partial result of algebraic manipulation:

More: How I would avoid most of the pesky algebra using the sympy library:

>>> R1, R2, R3, R4, R5, Vs1 = 1100, 2300, 1500, 550, 1600, 23
>>> from sympy import *
>>> var('I1,I2,I3,Vs2')
(I1, I2, I3, Vs2)
>>> eq1 = -Vs1 + R1*I1 + R2 * (I1-I2)
>>> eq1
3400*I1 - 2300*I2 - 23
>>> eq2 = R2*(I2-I1)+R3*I2+R4*(I2-I3)
>>> eq2
-2300*I1 + 4350*I2 - 550*I3
>>> eq3 = R4*(I3-I2)+R5*I3 + Vs2
>>> eq3
-550*I2 + 2150*I3 + Vs2
>>> from scipy import linalg
>>> import numpy as np
>>> for Vs2 in [10,15,20,25]:
... 	ls = np.array([[3400,-2300,0],[-2300,4350,-550],[0,-550,2150]])
... 	rs = np.array([23, 0, -Vs2])
... 	I = linalg.solve(ls, rs)
... 	Vs2, I
... 	
(10, array([ 0.01007914,  0.0048996 , -0.00339778]))
(15, array([ 0.00975305,  0.00441755, -0.00584667]))
(20, array([ 0.00942696,  0.0039355 , -0.00829557]))
(25, array([ 0.00910087,  0.00345346, -0.01074446]))

In order to solve a linear system of equations for unknown vector x=In which is classically written as Ax=b, you need to specify a coefficient matrix A and right hand side vector b to linalg.solve function. Based on your question, you just have to re-write in matrix form your three equations in terms of unknown currents to get A and b which was done with sympy but it is pretty overkill here imo. Here follows an easier to read solution with analytic A:

from scipy.linalg import lu_factor, lu_solve
import numpy as np
import matplotlib.pyplot as plt

# your data
R1 = 1100
R2 = 2300
R3 = 1500
R4 = 550
R5 = 1600
Vs1 = 23
# Vs2 range of interest as a list
Vs2_range = [10,15,20,25]

# construct A: the coefficient matrix of the left-hand side in terms of In = [I1, I2, I3]
A = np.array([[ R1+R2,       -R2,     0],
              [   -R2,  R2+R3+R4,   -R4],
              [     0,       -R4, R4+R5]])

# pre-compute pivoted LU decomposition of A to solve Ax=b (because only b is changing here)
A_LU = lu_factor(A)

# initialize results
res = np.empty((len(Vs2_range),3))

# loop over Vs2 values
for i,Vs2 in enumerate(Vs2_range):
    # construct b: the right hand side vector for each Vs2
    b = np.array([Vs1,0,-Vs2])
    # then solve the linear system Ax=b
    In = lu_solve(A_LU,b)
    # stock results as rows of the res array
    res[i,:] = In

# plot each current In (column of res) vs Vs2_range    
for i in range(3):    
    plt.plot(Vs2_range,res[:,i],'-+',label='I'+str(i+1))  
plt.xlabel('Vs2 [V]')
plt.ylabel('I [A]')
plt.legend()

which gives:

Hope this helps.