title: What type of (probably syntactic) mistake am I making when using a generic function on an array for plotting?
body: <p>I am trying to plot a geodesic on a 3D surface (tractrix) with Matlab. This worked for me in the past when I didn't need to parametrize the surface (see <a href="https://stackoverflow.com/a/60245213/4089351">here</a>). However, the tractrix called for parameterization, chain rule differentiation, and collection of <code>u,v,x,y</code> and <code>f(x,y)</code> values.</p>
  <p>After many mistakes I think that I'm getting the right values for <code>x = f1(u,v)</code> and <code>y = f2(u,v)</code> describing a spiral right at the base of the surface:</p>
  <p><a href="https://i.stack.imgur.com/5Prts.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/5Prts.png" alt="enter image description here" /></a></p>
  <p>What I can't understand is why the <code>z</code> value or height of the 3D plot of the curve is consistently zero, when I'm applying the same mathematical formula that allowed me to plot the surface in the first place, ie. <code>f  = @(x,y) a.* (y - tanh(y)) </code>.</p>
  <p>Here is the code, which runs without any errors on Octave. I'm typing a special note in upper case on the crucial calls. Also note that I have restricted the number of geodesic lines to 1 to decrease the execution time.</p>
  <pre><code>a = 0.3;
  
  u = 0:0.1:(2 * pi);
  v = 0:0.1:5;
  
  [X,Y] = meshgrid(u,v);
  % NOTE THAT THESE FORMULAS RESULT IN A SUCCESSFUL PLOT OF THE SURFACE:
  x = a.* cos(X) ./ cosh(Y);
  y = a.* sin(X) ./ cosh(Y);
  z = a.* (Y - tanh(Y));
  
  h = surf(x,y,z);
  zlim([0, 1.2]);
  set(h,'edgecolor','none')
  colormap summer
  hold on
  
  % THESE ARE THE GENERIC FUNCTIONS (f) WHICH DON'T SEEM TO WORK AT THE END: 
  f  = @(x,y) a.* (y - tanh(y)); 
  f1 = @(u,v) a.* cos(u) ./ cosh(v);
  f2 = @(u,v) a.* sin(u) ./ cosh(v);
  
  dfdu = @(u,v) ((f(f1(u,v)+eps, f2(u,v)) - f(f1(u,v) - eps, f2(u,v)))/(2 * eps) .*
                         (f1(u+eps,v)-f1(u-eps,v))/(2*eps) +
                         (f(f1(u,v), f2(u,v)+eps) - f(f1(u,v), f2(u,v)-eps))/(2 * eps) .*
                         (f2(u+eps,v)-f2(u-eps,v))/(2*eps));
  
  dfdv = @(u,v) ((f(f1(u,v)+eps, f2(u,v)) - f(f1(u,v) - eps, f2(u,v)))/(2 * eps) .*
                         (f1(u,v+eps)-f1(u,v-eps))/(2*eps) +
                         (f(f1(u,v), f2(u,v)+eps) - f(f1(u,v), f2(u,v)-eps))/(2 * eps) .*
                         (f2(u,v+eps)-f2(u,v-eps))/(2*eps));
  
   
  % Normal vector to the surface:
  N = @(u,v) [- dfdu(u,v), - dfdv(u,v), 1]; % Normal vec to surface @ any pt.
  
  % Some colors to draw the lines:
  C = {'k','r','g','y','m','c'};
  
  for s = 1:1     % No. of lines to be plotted.
    
  % Starting points:
  u0 = [0, u(length(u))];   
  v0 = [0, v(length(v))]; 
  du0 = 0.001; 
  dv0 = 0.001;  
     
  step_size = 0.00005; % Will determine the progression rate from pt to pt.
  eta =  step_size / sqrt(du0^2 + dv0^2); % Normalization.
  eps = 0.0001;          % Epsilon
  max_num_iter = 100000; % Number of dots in each line.
   
   
  % Semi-empty vectors to collect results:
  U = [[u0(s), u0(s) + eta*du0], zeros(1,max_num_iter - 2)]; 
  V = [[v0(s), v0(s) + eta*dv0], zeros(1,max_num_iter - 2)]; 
  
  
  for i = 2:(max_num_iter - 1)  % Creating the geodesic:
   
    ut = U(i); 
    vt = V(i);
    xt = f1(ut,vt);
    yt = f2(ut,vt);
    ft = f(xt,yt);
  
    utm1 = U(i - 1);  
    vtm1 = V(i - 1);
    xtm1 = f1(utm1,vtm1);
    ytm1 = f2(utm1,vtm1);
    ftm1 = f(xtm1,ytm1);
  
    usymp = ut + (ut - utm1);
    vsymp = vt + (vt - vtm1);
    xsymp = f1(usymp,vsymp);
    ysymp = f2(usymp,vsymp);
    fsymp = ft + (ft - ftm1);
  
    df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
    n = N(ut,vt);                % Normal vector at point t
    gamma = df * n(3);           % Scalar x change f x z value of N
  
    xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x &amp; y.
    ytp1 = ysymp - gamma * n(2);
  
    U(i + 1) = usymp - gamma * n(1);;
    V(i + 1) = vsymp - gamma * n(2);;
   
   end
  
   % THE PROBLEM! f(f1(U,V),f2(U,V)) below YIELDS ALL ZEROS!!! The expected values are between 0 and 1.2.
   P = [f1(U,V); f2(U,V); f(f1(U,V),f2(U,V))]; % Compiling results into a matrix.
  
      units = 35; % Determines speed (smaller, faster)
      packet = floor(size(P,2)/units);
      P = P(:,1: packet * units);
  
  
    for k = 1:packet:(packet * units)
          hold on
  
          plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))), 
              '.', 'MarkerSize', 5,'color',C{s})
          drawnow
    end
    end
  </code></pre>