How to setup MATLAB shape function solver

I am trying to find the shape function for a finite element analysis problem in MATLAB. The basic shape function looks like, with the boundary condition on the left. At Ni=N1, x=x1, y=y1 (the coordinate of node 1), the shape function would be equal to 1, when Ni=N1, but x=x2,y=y2 (the coordinate of node 2), the shape function would be 0, and so on and similar at Ni=N2,N3 and such. I use a symbolic 16*16 matrix named A to represent the 4 coefficients at the four shape functions.

This is my code to find the alpha, beta, gamma, and that last greek letter for the four equation from N1 to N4 (totally 16 unknown):

A = sym('A', [4 4]);
eqns = zeros(4,4);
coorx = sym('coorx', [1 4]);
coory = sym('coory', [1 4]);
for i=1:4
    for j=1:4
        if j==i
           eqns(i,j)=A(i,1)+A(i,2)*coorx(j)+A(i,3)*coory(j)+A(i,4)*coorx(j)*coory(j)==1;
        else
           eqns(i,j)=A(i,1)+A(i,2)*coorx(j)+A(i,3)*coory(j)+A(i,4)*coorx(j)*coory(j)==0;
        end
    end
end
tst=solve(eqns,A);
tst.A4_4

I tested out a lot of outputs, but the only value I got is 0 when I should get a function of sorts like coorx(j)-coory(j). Could you tell me why?

These are easy to write if you use parametric coordinates and shape functions. The usual order is to start at the lower left corner and proceed in a counterclockwise direction:

-1 <= a <= +1
-1 <= b <= +1

h1(a, b) = (1-a)*(1-b)/4.0
h2(a, b) = (1+a)*(1-b)/4.0
h3(a, b) = (1+a)*(1+b)/4.0
h4(a, b) = (1-a)*(1+b)/4.0

You use the Jacobean matrix to map from parametric to global coordinates.