numpy linear algebra solver

Let me get straight to the deal and describe the simple math problem I'm trying to solve in python:

I have 4 nxn stiffness matrices of coefficients, namely s11, s12, s21, s22, such that the system of algebraic equations I am trying to solve looks like the following:

[s11]*{u}+[s12]*{v} = {fx}
[s22]*{v}+[s21]*{u} = {fy}

where u and v - are nx1 components of unknown vector in the x and y directions. I.e. total vector will be g = sqrt(u^2+v^2). And fx and fy - are the nx1 x and y components of a given say "forces" vector.

So the task is relatively easy, but I'm stuck in implementing the gaussian elimination using the numpy in python. One of the first and most obvious methods I was trying to use, was rewriting everything in the following format:

s = [ [[s11], [s12]],
      [[s21], [s22]] ]
f = [ [fx],
      [fy] ]

Such that s is 2nx2n and f is 2nx1. The way of assembling them is the following:

s1 = np.concatenate((s11, s12), axis=1)
s2 = np.concatenate((s21, s22), axis=1)
s = np.concatenate((s1, s2), axis=0)
f = np.concatenate((fx, fy), axis=0)

But then, when I try to solve the system with

u = np.linalg.solve(s, f)

I expect to get 2nx1 vector of unknowns u (where the first set of n elements are the x components and second set is y), but I get the following error:

Traceback (most recent call last):
File "/home/art/PycharmProjects/Cavern/main.py", line 204, in <module>
u = np.linalg.solve(s, f)
File "<__array_function__ internals>", line 5, in solve
File "/home/art/.pyenv/versions/test380/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 399, in solve
r = gufunc(a, b, signature=signature, extobj=extobj)
File "/home/art/.pyenv/versions/test380/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 97, in _raise_linalgerror_singular
raise LinAlgError("Singular matrix")
numpy.linalg.LinAlgError: Singular matrix

Then I have tried another controversial way, which can be summarized as follows:

s1 = s11 + s21
s2 = s12 + s22
s = np.concatenate((s1, s2), axis=1)
f = fx + fy

Such that s is now nx2n, f is nx1. Then again, I'm trying to use the gaussian elimination as follows:

u = np.linalg.solve(s, f)

Again I expect to get 2nx1 vector of unknowns, but I get the error which now is the following:

Traceback (most recent call last):
File "/home/art/PycharmProjects/Cavern/main.py", line 204, in <module>
u = np.linalg.solve(s, f)
File "<__array_function__ internals>", line 5, in solve
File "/home/art/.pyenv/versions/test380/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 386, in solve
_assert_stacked_square(a)
File "/home/art/.pyenv/versions/test380/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 213, in _assert_stacked_square
raise LinAlgError('Last 2 dimensions of the array must be square')
numpy.linalg.LinAlgError: Last 2 dimensions of the array must be square

Therefore I keen to seek for your help. Any assistance/hints would be much appreciated. And yes, mb it's not important to mention, but the whole thing is about transferring my code from Matlab to python. In the Matlab the procedure is the following:

S1 = [S11, S12]; S2 = [S21, S22];
U1 = S1 \ fx; U2 = S2 \ fy;
U = U1 + U2;
u(:,1) = U(1:nnodes);
v(:,1) = U(nnodes+1:end);

Backslash stands for the gaussian elimination here.

Why use np.concatenate if the s and f coefficients are scalars? You can directly make them into numpy arrays (or any ordered sequence) and use numpy.linalg.solve() on them.

s1 = [s11, s12]
s2 = [s21, s22]
s = np.column_stack((s1, s2))
f = [s11, s12]

u = np.linalg.solve(s, f)

Note in this case, the 1D list is the same as having a nx1 (or 1xn) vector. The directionality of the vector does not matter for the np.linalg.solve() function.