Entries doesn't have to be 1 to 16, it can be anything. As long as all numbers are different.

Here's a brute force solution.  To reduce the solution space, I solved for a pan-magic square; in addition to the usual constraints, the sums of all the broken diagonals must also equal the magic number.

PanMagic := proc(L,n::posint)
local i,j,M,rng;
    rng := 0..n-1;
    M := (n^3+n)/2;
    [seq(s=M, s in [NULL
                    , seq(add(L[i,j], j=rng), i=rng)
                    , seq(add(L[i,j], i=rng), j=rng)
                    , seq(add(L[i,modp(i+j,n)],i=rng), j=rng)
                    , seq(add(L[i,modp(j-i-1,n)],i=rng), j=rng)
                   ])];
end proc:

to 1 do
    n := 4:
    eqs := PanMagic(L,n);
    sol := solve(eqs);
    vals := map(rhs,sol);
    vars := convert(indets(vals),'list');
end do:

all := {seq(1..n^2)}:

P := combinat:-permute(convert(all,'list'),nops(vars)):

found := false:
for p in P do
    eqs := Equate(vars,p);
    if subs(eqs, vals) = all then
        found := true;
        break;
    end if;
end do:

if found=true then
    M := subs(sol,eqs, Matrix(n,n,(i,j) -> L[i-1,j-1]));
end if;

                               [ 1     8    10    15]
                               [                    ]
                               [12    13     3     6]
                          M := [                    ]
                               [ 7     2    16     9]
                               [                    ]
                               [14    11     5     4]

You seem to be insisting on doing things the hard way.  Moreover, you're using the old and deprecated linalg package rather than LinearAlgebra.  But anyway:
 

> PL := convert(P, list);
   vars:= map(parse, [$"a".."p"]);

Get some arbitrary values for v and _t[1] to _t[7], and get the variable values:

> vals:= subs(map(t -> (t = rand(0..100)()), indets(PL)),zip(`=`,vars,PL));

> subs(vals, Matrix(4,4,vars));