Constructs an interpolating spline.
| xy | It is a matrix the first column of which contains the x values, the second contains the corresponding y values. The 'x' coordinates are assumed to be increasing. | 
| cond | (Optional) Additional (linear) constraints. For details see description. | 
| return value | Spline matching the spline and additional constraints. | 
The order of the spline will be the number of additional constraints plus one. In particular, (explicit) no constraints will imply a piece-wise linear function. However, the default is cubic interpolation with zero second derivatives at the two end points.
  xy = [-100,4; 8,93; 42,-88; 444, 827];
  xunit = 1J;
  yunit = 1V;
  xmin = min(xy[...;0]);
  xmax = max(xy[...;0]);
  eps = (xmin-xmax)/200;
  xmin -= 2*eps;
  xmax += 2*eps;
  x = (xmin..eps..xmax)*xunit;
  print("Default constraints").title()
  p = splineInterp(xy * [xunit,0;0,yunit])
  y = (@(x) splineEval(p,x))(#x);
  [! xy[...;0]*xunit, (@(x) splineEval(p,x))(#(xy[...;0]*xunit))' !]
  plot(x,y)
  print("Piece-wise linear").title()
  p = splineInterp(xy * [xunit,0;0,yunit], [])
  y = (@(x) splineEval(p,x))(#x);
  plot(x,y)
  print("Random constaints").title()
  cond = [
      // It will take 4 at 128 (but 128 is not a node for the spline)
      (@(s) splineEval(s,128J)), 4V;
      // The second derivative at 228 will be also 4 (without being a node)
      (@(s) splineEval(splineDiff(s,2),228J)*1J^2), 4V;
      // ...and some combination of the above two.
      (@(s) splineEval(splineDiff(s,3), 16J)*1J^3 - 2*splineEval(s,128J)), 42V];
  p = splineInterp(xy * [xunit,0;0,yunit], cond)
  y = (@(x) splineEval(p,x))(#x);
  plot(x,y)