from gekko import GEKKO

m = GEKKO()

ni = 3  # number of rows
nj = 2  # number of columns

# best method: use m.Array function
x = m.Array(m.Var,(ni,nj))
m.Equations([x[i][j]==i*j+1 for i in range(ni) for j in range(nj)])

# another way: list comprehensions
y = [[m.Var() for j in range(nj)] for i in range(ni)]
for i in range(ni):
     for j in range(nj):
         m.Equation(x[i][j]**2==y[i][j])

# summation
z = m.Var()
m.Equation(z==sum([sum([x[i][j] for i in range(ni)]) for j in range(nj)]))

m.solve()

print('x:')
print(x)
print('y=x**2:')
print(y)
print('z')
print(z.value)