Algorithm/Code for finding surface shape of an image?
The 3x3 linear system appears to be that for photometric stereo, where the $i$-th row is $$l_i\cdot n = b_i,$$
with


$l_i$ is the 3D light source vector, e.g., $l_1$ the first row of the matrix, (0.2425, 0, -0.9701)
$n = (x_1, x_2, x_3)$ is the unknown surface normal vector, or "shape"
$b_i$ is the brightness observed when the surface was illuminated with $l_i$


Solving this 3x3 linear system gives the normal $n$ for a single surface point ("middle" according to "Example 1").  (You may wish to divide $n$ by its length $\sqrt{x_1^2+x_2^2+x_3^2}$ so that it is unit length.)  

Usually, by surface shape in this context of photometric stereo one means the height function $z(x,y)$, where we associate a height to each pixel $(x,y)$.  To obtain this, we first need to solve a related 3x3 system at each pixel location $(x,y)$ (not just the "middle"), so that we have a surface normal $n(x,y)$ everywhere.  Then we solve for that function $z$ whose normal function is that $n$ we obtained via many 3x3 solutions above.  Finding $z$ requires we impose integrability on vector field $n$ (because modeling error and noise may mean there exists no $z$ having normals $n$).  However, I am digressing from your original "Example 1", which seeks a single surface normal.

See also: http://pages.cs.wisc.edu/~lizhang/courses/cs766-2008f/syllabus/10-09-shading/shading.pdf and https://en.wikipedia.org/wiki/Photometric_stereo

Note that the above assumes diffuse reflection (no surface glossiness), surface smoothness, and no shadows nor occlusion.  Light sources must be distant if we keep them constant over image.