|   | 1 |  initial version  | 
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
 Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
 
                
                Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.