Here is a way to build a linear system as in the question.

First, define a function `a`

to create symbolic variables $a_{i,j}$.

Here we represent $a_{i, j}$ as `a_i_j`

, and we use `m`

for the
minus sign, so `m1`

stands for $-1$.

```
def a(t):
return SR.var(('a' + '_{}'*len(t)).format(*[str(i).replace('-','m') for i in t]))
```

Test it:

```
sage: a((1, 1)), a((12, 437)), a((-1, 2)), a((-1, 2, 8)), a((2, 4, -1, 3))
(a_1_1, a_12_437, a_m1_2, a_m1_2_8, a_2_4_m1_3)
```

Then, we need a function to build the linear system in the question.

```
def thesystem(p, q):
s = []
for i in range((p-1)//2 + 1):
for j in range(q-2):
s.append(a((i, j)) - 3*a((i-1, j)) - 4*a((i, j-1)) + 10*a((i-1, j-1)) == 0)
for j in range(q-1):
s.append(a((-1, j)) + a(((p-1)//2 - 1, j)) == 0)
for i in range((p-1)//2 + 1):
s.append(a((i, -1)) - a((i, q-2)) == 0)
return s
```

Test it with $(p, q) = (5, 4)$:

```
sage: s_5_4 = thesystem(5, 4)
sage: s_5_4
[a_0_0 - 4*a_0_m1 - 3*a_m1_0 + 10*a_m1_m1 == 0,
-4*a_0_0 + a_0_1 + 10*a_m1_0 - 3*a_m1_1 == 0,
-3*a_0_0 + 10*a_0_m1 + a_1_0 - 4*a_1_m1 == 0,
10*a_0_0 - 3*a_0_1 - 4*a_1_0 + a_1_1 == 0,
-3*a_1_0 + 10*a_1_m1 + a_2_0 - 4*a_2_m1 == 0,
10*a_1_0 - 3*a_1_1 - 4*a_2_0 + a_2_1 == 0,
a_1_0 + a_m1_0 == 0,
a_1_1 + a_m1_1 == 0,
a_1_2 + a_m1_2 == 0,
-a_0_2 + a_0_m1 == 0,
-a_1_2 + a_1_m1 == 0,
-a_2_2 + a_2_m1 == 0]
```

Now a function to extract the variables from the linear system:

```
def sysvar(s):
v = set()
for e in s:
v.update(e.variables())
return sorted(v, key=str)
```

Try it:

```
sage: v_5_4 = sysvar(s_5_4)
sage: v_5_4
[a_0_0,
a_0_1,
a_0_2,
a_0_m1,
a_1_0,
a_1_1,
a_1_2,
a_1_m1,
a_2_0,
a_2_1,
a_2_2,
a_2_m1,
a_m1_0,
a_m1_1,
a_m1_2,
a_m1_m1]
```

Now we can solve the linear system:

```
sage: solve(s_5_4, v_5_4)
[[a_0_0 == r8,
a_0_1 == -1/5*r5 + 17/5*r8,
a_0_2 == -1/10*r5 + 2/5*r7 + 3/10*r8,
a_0_m1 == -1/10*r5 + 2/5*r7 + 3/10*r8,
a_1_0 == r5,
a_1_1 == 17/5*r5 + 1/5*r8,
a_1_2 == r7,
a_1_m1 == r7,
a_2_0 == -1/20*r5 + 1/4*r6 - 3/20*r8,
a_2_1 == r6,
a_2_2 == -61/80*r5 + 1/16*r6 + 5/2*r7 - 3/80*r8,
a_2_m1 == -61/80*r5 + 1/16*r6 + 5/2*r7 - 3/80*r8,
a_m1_0 == -r5,
a_m1_1 == -17/5*r5 - 1/5*r8,
a_m1_2 == -r7,
a_m1_m1 == -17/50*r5 + 4/25*r7 + 1/50*r8]]
```

The solution involves some real parameters `rj`

whose numbering
may vary depending on the history of your Sage session.

Alternatively, we could build a matrix for the linear system.

For this, we define a function, inspired by
this answer
to this question,
but taking advantage of the fact that there are no constant terms here.

```
def matrix_from_system(s, v):
return matrix([[(e.lhs() - e.rhs()).coefficient(a) for a in v] for e in s])
```

The following matrix is obtained by applying this function to our example:

```
sage: matrix_from_system(s_5_4, v_5_4)
[ 1 0 0 -4 0 0 0 0 0 0 0 0 -3 0 0 10]
[-4 1 0 0 0 0 0 0 0 0 0 0 10 -3 0 0]
[-3 0 0 10 1 0 0 -4 0 0 0 0 0 0 0 0]
[10 -3 0 0 -4 1 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 -3 0 0 10 1 0 0 -4 0 0 0 0]
[ 0 0 0 0 10 -3 0 0 -4 1 0 0 0 0 0 0]
[ 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0]
[ 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0]
[ 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0]
[ 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0]
```

and we can then use linear algebra.

For more on linear algebra, the approach using matrices and vectors,
and the approach of solving linear systems, see

Is $p$ assumed to be an odd integer and $\ge 1$, so that $\frac{p-1}{2}$ is a non-negative integer?

Is $q$ assumed to be an integer and $\ge 2$, so that $q-2$ is a non-negative integer?