This looks pretty good! I've left some comments based on a quick read (did not check the math yet).
I'd like to see a test with random values for x and y, and some testing with non-float input. Also handling of input with nans should be fixed in tests (there will be user requests :-).

If this is supposed to be beginner-friendly, some more commenting of the construction would be nice. A simple example in the docstring would be nice, too.

As far as the math is concerned, it is worth spending a little bit of effort to check the numerical stability of the specific implementation.

Very useful notes, Evgeni. I addressed most of the things.

What about NaN's handling?

As example I'm thinking about the plot of derivatives 0, 1, 2, 3 to show that the continuity property breaks only for the third derivative.

NaN handling: I'd suggest to 1) add a test that data with nans raises, and 2) add a note to the docstring. Expect a bunch of pandas users who expect some sort of nan-means-missing-value intelligence :-).

As example I'm thinking about the plot of derivatives 0, 1, 2, 3 to show that the continuity property breaks only for the third derivative.

+1

I came up with a way to check continuity of the solution, using 1st order Taylor approximation to compare close values near break points. This way I could sensibly check solution properties with x and y randomly generated.

Regarding the first test: one can obtain the same values in Octave using the code provided.

An example is here. More or less it is ready, so I added [MRG] tag.

@@            master   #5653   diff @@
======================================
  Files          235     235       
  Stmts        43282   43347    +65
  Branches      8155    8164     +9
  Methods          0       0       
======================================
+ Hit          33664   33723    +59
- Partial       2591    2594     +3
- Missed        7027    7030     +3

Review entire Coverage Diff as of 9ac3b61

Powered by Codecov. Updated on successful CI builds.

This definitely needs an entry in the release notes.

I wrote a note to release notes.

Fixed problems after rebase.

@ev-br yeah, ok.

I reverted changes, look forward to see a new PR from Matthieu

CI failure is due to the bot complaining about whitespace.

$ pep8 scipy benchmarks/benchmarks
benchmarks/benchmarks/signal_filtering.py:30:5: E303 too many blank lines (2)
1       E303 too many blank lines (2)

BTW, gh-5734 in part addresses concerns of @juliantaylor and @argriffing on what is the value added by CubicSpline as compared by UnivariateSpline(x, y, s=0). Flexibility in specifying boundary conditions is definitely a useful feature, and adding it to fitpack-based splines is awkward.

I plan to merge this soonish, unless further comments.

needs a rebase now

Rebased.

Merged.
This is going to be one of the highlights of the release, I'd think. Thanks Nikolay, and keep them coming!