One way would be to use fzero. If the parametric equations for your line are L(t)=L0+t*D, the t that gives the intersection could be obtained as,
t_intersection = fzero(@(t) objective(t,L0,D,fitobj), t0)
function difference=objective(t,L0,D,fitobj)
L=L0+t*D;
x=L(1); y=L(2); z=L(3);
difference = z-fitobj(x,y);
end
