Ti - basic homography bilinear interpolation problem

I'm trying to make a program that can calculate bilinear interpolation with homography; which means I have four points (xn,yn) n = 1, 2, 3, 4, and corresponding f(xn, yn) and interpolate f(x,y) from the given four sets.

The basic concept is to get a homography matrix that moves four given (xn, yn) to the unit square, and to apply bilinear interpolation to the unit square.

https://math.stackexchange.com/questions/494238/how-to-compute-homography-matrix-h-from-corresponding-points-2d-2d-planar-homog (see Vishnu Pardeesh's answer)

https://en.wikipedia.org/wiki/Bilinear_interpolation

https://www.cs.ubc.ca/grads/resources/thesis/May09/Dubrofsky_Elan.pdf

And below is my Ti-basic code(including debugging codes)

Define LibPub blin()=
Prgm
:Request "x1",x1
:Request "y1",y1
:Request "f1",f1
:Request "x2",x2
:Request "y2",y2
:Request "f2",f2
:Request "x3",x3
:Request "y3",y3
:Request "f3",f3
:Request "x4",x4
:Request "y4",y4
:Request "f4",f4
:Request "x",x
:Request "y",y
:Local a
:a:=[[−x1,−y1,−1,0,0,0,0,0,0][0,0,0,−x1,−y1,−1,0,0,0][−x2,−y2,−1,0,0,0,x2,y2,1][0,0,0,−x2,−y2,−1,0,0,0][−x3,−y3,−1,0,0,0,x3,y3,1][0,0,0,−x3,−y3,−1,0,0,0][−x4,−y4,−1,0,0,0,x4,y4,1][0,0,0,−x4,−y4,−1,x4,y4,1][0,0,0,0,0,0,0,0,1]]
:Local b
:b:=[[0][0][0][0][0][0][0][0][1]]
:Local h
:h2:=(a*a)^(−1)*a*b
:Disp "h2 = ",h2
:Disp "b=a*h=",a*h2
:h:=[[h2[1,1],h2[2,1],h2[3,1]][h2[4,1],h2[5,1],h2[6,1]][h2[7,1],h2[8,1],h2[9,1]]]
:Disp "Homography h =",h
:Local q,q1,q2,q3,q4
:q:=h*[[x][y][1]]
:q1:=h*[[x1][y1][1]]
:q2:=h*[[x2][y2][1]]
:q3:=h*[[x3][y3][1]]
:q4:=h*[[x4][y4][1]]
:q:=((q)/(q[3,1]))
:q1:=((q1)/(q1[3,1]))
:q2:=((q2)/(q2[3,1]))
:q3:=((q3)/(q3[3,1]))
:q4:=((q4)/(q4[3,1]))
:Disp q,q1,q2,q3,q4
:Local f
:f:=[1-q[1,1],q[1,1]]*[[f1,f3][f2,f4]]*[[1-q[2,1]][q[2,1]]]
:Disp "f=",f
:Return f
:EndPrgm

And the testing result is like below.

x1 140
h2 =  [[−0.007144][0.182382][1.][7.70819−14][1.08176−12][−9.54405−12][−0.007144][0.182382][1.]]
b=a*h= [[1.77−11][−1.24859−12][1.57−11][−1.79793−12][1.86−11][−2.01942−12][1.78−11][−1.44315−12][1.]]
Homography h = [[−0.007144,0.182382,1.][7.70819−14,1.08176−12,−9.54405−12][−0.007144,0.182382,1.]]
[[1.][−3.2999−10][1.]] [[17.7][−1.24859][1.]] [[1.][1.94125−11][1.]] [[1.][−2.82671−11][1.]] [[−16.8][2.44315][1.]]
f= [2550.]

Since some fonts are breaking, I append screenshot here.

It seems the homography matrix does not move the given points to a unit square. How can I fix the code?

It turns out that the order of the last three columns of the fifth and sixth row of the A has been changed.

thank you.