Alternative way is expressing these two line by y-y0=k(x-x0) . And solve the linear system as you did. this code didn't handle degenerate configurations.

proc iml;
start IntersectSegsSimple(p1, p2, q1, q2);
k1=(p1[2]-p2[2])/(p1[1]-p2[1]);
k2=(q1[2]-q2[2])/(q1[1]-q2[1]);
A = (k1//k2)||{-1,-1};
b = (k1#p1[1]-p1[2])//(k2#q1[1]-q1[2]);
x = solve(A, b);
return (x);
finish;
/* Test 1: intersection at (0.95, 1.25) */
p1 = {1.8 2.1}; p2 = {0.8 1.1};
q1 = {1 1.25}; q2 = {0 1.25};
z = IntersectSegsSimple(p1,p2,q1,q2);
print z;
quit;