How do I extract the intersection line between a plane and a surface?

I have 2 sets of 3D points through which I fitted a plane and a surface (see attached figure). I tried to extract the line of intersection using symbolic expressions as described in: https://uk.mathworks.com/matlabcentral/answers/301302-intersection-of-two-sfit-planes However, I am unable to visualize the intersection line as I get the error “Error using inlineeval (line 14) Error in inline expression… “. More importantly, I’m stuck as to how to get the coordinates of x,y,z points along the intersection line.

clc; clear; close all;
       
Xs = [18.5, 18, 17.5
      18.5, 18, 16.5
      19.5, 18, 17.5
      19.5, 18, 16.5
      17.5, 18, 18.5
      16.5, 18, 18.5];
 
Xn = [18.5, 18.5, 18.0
      18.5, 19.5, 18.0
      19.5, 18.5, 18.0
      19.5, 19.0, 17.5
      19.5, 19.5, 17.0
      18.0, 18.5, 18.5
      18.0, 19.5, 18.5];
  
% fit plane through Xs
[sfp,gofp,obj] = fit([Xs(:,1), Xs(:,3)],Xs(:,2), 'poly11');
plot(sfp,[Xs(:,1), Xs(:,3)],Xs(:,2))
hold on;
 
% fit surface through Xn
[sfs,gofs,obj] = fit([Xn(:,1), Xn(:,3)],Xn(:,2), 'poly22');
plot(sfs,[Xn(:,1), Xn(:,3)],Xn(:,2))
 
% obtain coefficients for yplane and ysurface
g = sfp.p10;
h = sfp.p01;
k = sfp.p00;

a = sfs.p20;
b = sfs.p02;
c = sfs.p10;
d = sfs.p01;
e = sfs.p11;
f = sfs.p00;
 
syms x z
yp = g*x + h*z + k;
ys = a*(x^2) + b*(z^2) + c*x + d*z + e*x*z + f;
 
% Solve for z in the expression yp==ys. This gives you z(x). 
z_of_x = solve(yp==ys, z);
 
% Noow substitute the value of z(x) in yp to obtain an expression for y(x)
y_of_x = simplify(subs(yp, z, z_of_x));
 
% Plot the result
figure
h1 = ezsurf(yp,[-20 20]);
set(h1,'facecolor','r','facealpha',0.5);
hold  all
h2 = ezsurf(ys,[-20 20]);
set(h2,'facecolor','b','facealpha',0.5);
hp = ezplot3(x, z_of_x, y_of_x, [-2 2]);
set(hp,'LineWidth',2,'Color','g');
ylim([-2 2]);

[madhan ravi]

you can select the code and press the code button so that it's easy to read