1. The local coordinates should go from -1..1 (not 0..1).
2. Switch order for the local coordinates 1+/-s/t (not s/t+/-1).
3. Use the elemetwise product operator ".*" between s/t (s and t are matrices and using "*" will result in matrix multiplication which is not correct).
[s,t] = meshgrid(-1:0.1:1,0:0.1:1);
N1 = 1/4*(1-t).*(1-s);
subplot(2,2,1)
surf(s,t,N1)
N2 = @(s,t)1/4*(1-t).*(1+s);
subplot(2,2,2)
surf(s,t,N2(s,t))
N3 = @(s,t)1/4*(1+t).*(1+s);
subplot(2,2,3)
surf(s,t,N3(s,t))
N4 = @(s,t)1/4*(1+t).*(1-s);
subplot(2,2,4)
surf(s,t,N4(s,t))
