sys := ([seq])(seq(eq[i, j], j = 0 .. N), i = 1 .. N):

nops(sys);

vars:=indets(sys) minus {k}:

nn := Matrix(N+1, N+1,(i, j)-> u[i-1, j-1]):

##

p:=proc(kk) local u_res,A;

  u_res:=solve(eval(sys,k=kk),vars);

  A:=eval(nn,u_res);

  plots:-matrixplot(A)

end proc;

## Testing p for k=0.001:

p(0.001);

## Animating the plot for k=0.0001..0.001:

plots:-animate(p,[k],k=0.0001..0.001);