Finding a point on a surface given an arc length and direction/angle in Matlab

I need to find a point on a surface (given an angle relative to a starting point) where the arc length is a given value using Matlab.

Say I have a high order surface where z=f(x,y) which was been fitted from sampling points using Matlabs fit function. If I have a starting point, say a = (x_0, y_0, f(x_0,y_0)) and want to know the coordinate of a point along that surface at a user defined angle theta away in the xy plane so that the distance covered over the surface is a given value, e.g. 10mm.

I assume that what I need to do is solve this equation for the value of b, given we know a, s and function defining surface. But I'm unsure how to write this in Matlab. I'm assuming I need to use the solve function in Matlab.

Any help on how to write this in Matlab and in the most efficent form would be greatly appreciated!

here it is an example assuming dx=1,dy=1, incorporating arbitrary x and y step sizes should not be hard

% //I am assuming here that you know how to get your Z
z=peaks(60); 
% //start point
spoint=[30,30];
% //user given angle
angle=pi/4;
% // distance you want
distance=10;
%// this is the furthes the poitn can be
endpoint=[spoint(1)+distance*cos(angle) spoint(2)+distance*sin(angle)]; 

%// we will need to discretize, so choose your "accuracy"
npoints=100;
%//compute the path integral over the line defined by startpoitn and endpoint
[cx,cy,cz]=improfile(z,[spoint(1) endpoint(1)],[spoint(2) endpoint(2)],npoints);

% // this computes distances between adjacent points and then computes the cumulative sum
dcx=diff(cx);
dcy=diff(cy);
dcz=diff(cz);

totaldist=cumsum(sqrt(dcx.^2+dcy.^2+dcz.^2));
%// here it is! the last index before it gets to the desired distance
ind=find(totaldist<distance,1,'last');

Visualization code

imagesc(z);axis xy;colormap gray
hold on;
plot(spoint(1),spoint(2),'r.','markersize',10)
plot(endpoint(1),endpoint(2),'r*','markersize',5)
plot([spoint(1) endpoint(1)],[spoint(2) endpoint(2)],'b')
plot(cx(ind),cx(ind),'g*','markersize',10)