interp3() is supposed to do interpolation for volumetric data points (Xi, Yi, Zi, Si), i=1~N. It will generate a function S=f(X, Y, Z). This is the reason why it needs the additional input 's' as noted in your post.
When you encounter an interpolation problem/algorithm, the first thing you need to figure out is what kind of function you are looking for. Namely, whether the function is univariate (y=f(x)), bivariate (z=f(x,y)) or multivariate ( s=f(x, y, z, ....) ). For your particular problem where you want to interpolate a series of 3D points using a spline, basically it is a univariate interpolation problem. However, since a space curve cannot be represented as y=f(x), the spline function will be represented in parametric form as S(t)=(x(t), y(t), z(t)).
There are many ways for doing spline interpolation thru 3D data points. Among them, the two algorithms that are very easy to implement are Catmull Rom spline and Overhauser spline. Both are cubic spline and differ only in how the first derivatives at the data points are estimated. You can google them to find out the details.
Best Answer
You can't calculate the gradient of a "kink" (as you so eloquently put it). If you really need a gradient at such a point (x), I'd just average the gradient at (x-d) and (x+d) where d is a small enough delta. It's as mathematically valid as any other single answer you're likely to get.
For example, the function:
will produce:
where there is no gradient at the origin (0,0). However, averaging the gradients at -0.0001 (gradient = -1) and +0.0001 (gradient = +1) will give you a gradient of zero (flat line).
This should give a half-decent answer for other equations that produce non-symmetrical gradients at (x-d) and (x+d) as well.
What I would do, since it's licensed under MIT, is to modify the source to allow an option for the Bezierspline to use that +/- delta method to calculate gradients at the non-continuous points. Maybe even push back the source changes to the developers if you think it's a worthwhile addition.