You have an 2D mesh of height values. I'm going to assume you create the faces of the surface for the $i$,$j$th square using two triangles of $\left\langle\left(i,j\right),\left(i+1,j\right),\left(i+1,j+1\right)\right\rangle$ and $\left\langle\left(i,j\right),\left(i,j+1\right),\left(i+1,j+1\right)\right\rangle$.
Your vertical cut plain can be defined as a 2D line called $l$. $l$ can be represented by using 2 known points on the edge of the mesh or by using an 2D line equation $i = mj+c$, $l=\left\langle{m,c}\right\rangle$.
The intersection points are the intersection with any of the edges of triangles along that 2D line $l$.

All the points heights are worked out by using linear interpolation between the heights given by the two array coordinated connected by the triangle edge being intersected.
Distance along the cut plane is given by Pythagorean theorem.
#!/usr/bin/env python
import math
from random import choice, uniform
mesh = [[10,5,8],[12,7,5],[4,7,3],[4,6,4],[8,3,1]]
print mesh
iMax = len(mesh)-1
jMax = len(mesh[0])-1
# pick random vertical cut plane
if choice([True,False]):
p1 = (uniform(0,iMax), choice([0,jMax]))
else:
p1 = (choice([0,iMax]), uniform(0,jMax))
while True:
if choice([True,False]):
p2 = (uniform(0,iMax), choice([0,jMax]))
else:
p2 = (choice([0,iMax]), uniform(0,jMax))
if p1[0] != p2[0] and p1[1] != p2[1]:
break
print "P1",p1
print "P2",p2
if p1[1] == p2[1]:
m = float('Inf') # Not handling this, needs fix
c = None
else:
m = float(p2[0]-p1[0])/(p2[1]-p1[1])
c = p1[0]-m*p1[1]
# i=mj+c formular
print "i={0}j+{1}".format(m,c)
def floor(a):
return int(math.floor(a))
def ceil(a):
return int(math.ceil(a))
iL = floor(min(p1[0], p2[0]))
iU = ceil(max(p1[0], p2[0]))
jL = floor(min(p1[1], p2[1]))
jU = ceil(max(p1[1], p2[1]))
print "points to plot unordered, d is distance from p1 (pythagoras), h is height as given by linear interpolation"
def linearInterpolation(h0,h1,r):
return (h1-h0)*r+h0
def pythag(a,b):
return math.sqrt(a**2 + b**2)
for i in xrange(iL, iU+1):
j = (i-c)/m
if 0 <= j <= jMax:
j1, j2 = floor(j), ceil(j)
h = linearInterpolation(mesh[i][j1], mesh[i][j2], j - j1)
d = pythag(p1[0]-i, p1[1]-j)
print (d,h)
for j in xrange(jL, jU+1):
i = m*j+c
if 0 <= i <= iMax:
i1, i2 = floor(i), ceil(i)
h = linearInterpolation(mesh[i1][j], mesh[i2][j], i - i1)
d = pythag(p1[0]-i, p1[1]-j)
print (d,h)
for v in xrange(-iU-1, jU):
i = (c - v)/(1.0-m)
j = m*i+c
if 0 <= i <= iMax and 0 <= j <= jMax:
i1, i2 = floor(i), ceil(i)
j1, j2 = floor(j), ceil(j)
r = pythag(i - i1, j - j1) / math.sqrt(2)
h = linearInterpolation(mesh[i1][j1], mesh[i2][j2], r)
d = pythag(p1[0]-i, p1[1]-j)
print (d,h)