MATLAB: Custom cubic spline function giving undesirable jumps in the derivative and double derivative. How to solve the issue

cubic splinescurve fittingCurve Fitting Toolbox

In the figure below you see the following:
C1 C2
R1 data1 data2
R2 d d
R3 dd dd
d: derivative
dd double derivative
red: original data
blue: cubic spline fit
It can be seen that the derivative and double derivative have impulsive jumps which should not be the case theoretically for splines, as it ensures C0, C1 and C2 continuity.
Below you can find the custom code which solves the following system of linear eqyation: Ma=y; a=M\y. This can also be seen in the attachment of splines here.
What is the reason for such behaviour and how can I eliminate that?
Custom Code:
function output = cubic_spline(X,Y,dX)%datax, datay, endslopes
n = length(X);
output.form = 'pp';
output.breaks = X';
output.coefs = zeros(n-1,4);
output.pieces = n-1;
output.order = 4;
output.dim = 1;
% Form: Ma=y, a:coefficients
M = zeros((n-1)*4,(n-1)*4);
y = zeros((n-1)*4,1);
c0 = [1 1 1 1];
c1 = [0 1 2 3];
c2 = [0 0 2 6];
%Starting slope end condition
M(1,1) = 0; M(1,2) = X(1,1)^0; M(1,3) =2*X(1,1)^1; M(1,4) = 3*X(1,1)^2;
y(1,1) = dX(1);
for i=2:1+(n-2) %For n-2 segments
    %Starting continuity condition

    M((i-2)*4+2,((i-2)*4)+1) = X(i-1,1)^0; M((i-2)*4+2,((i-2)*4)+2) = X(i-1,1)^1; M((i-2)*4+2,((i-2)*4)+3) = X(i-1,1)^2; M((i-2)*4+2,((i-2)*4)+4) =  X(i-1,1)^3;
    y((i-2)*4+2,1) = Y(i-1,1);
    %C0 continuity conditions

    M((i-2)*4+3,((i-2)*4)+1) = X(i,1)^0; M((i-2)*4+3,((i-2)*4)+2) = X(i,1)^1; M((i-2)*4+3,((i-2)*4)+3) = X(i,1)^2; M((i-2)*4+3,((i-2)*4)+4) = X(i,1)^3;
    y((i-2)*4+3,1) = Y(i,1);
    %C1 continuity conditions
    M((i-2)*4+4,((i-2)*4)+1) = 0; M((i-2)*4+4,((i-2)*4)+2) = X(i,1)^0; M((i-2)*4+4,((i-2)*4)+3) = 2*X(i,1)^1; M((i-2)*4+4,((i-2)*4)+4) = 3*X(i,1)^2; M((i-2)*4+4,((i-2)*4)+5) = -0; M((i-2)*4+4,((i-2)*4)+6) = -X(i,1)^0; M((i-2)*4+4,((i-2)*4)+7) = -2*X(i,1)^1; M((i-2)*4+4,((i-2)*4)+8) = -3*X(i,1)^2;
    y((i-2)*4+4,1) = 0;
    %C2 continuity conditions
    M((i-2)*4+5,((i-2)*4)+1) = 0; M((i-2)*4+5,((i-2)*4)+2) = 0; M((i-2)*4+5,((i-2)*4)+3) = 2*X(i,1)^0; M((i-2)*4+5,((i-2)*4)+4) = 6*X(i,1)^1; M((i-2)*4+5,((i-2)*4)+5) = -0; M((i-2)*4+5,((i-2)*4)+6) = -0; M((i-2)*4+5,((i-2)*4)+7) = -2*X(i,1)^0; M((i-2)*4+5,((i-2)*4)+8) = -6*X(i,1)^1;
    y((i-2)*4+5,1) = 0; 
end
%Last segment continuity condition
i = i+1;
    %Starting continuity condition
    M((i-2)*4+2,((i-2)*4)+1) = X(i-1,1)^0; M((i-2)*4+2,((i-2)*4)+2) = X(i-1,1)^1; M((i-2)*4+2,((i-2)*4)+3) = X(i-1,1)^2; M((i-2)*4+2,((i-2)*4)+4) =  X(i-1,1)^3;
    y((i-2)*4+2,1) = Y(i-1,1);
    %C0 continuity conditions
    M((i-2)*4+3,((i-2)*4)+1) = X(i,1)^0; M((i-2)*4+3,((i-2)*4)+2) = X(i,1)^1; M((i-2)*4+3,((i-2)*4)+3) = X(i,1)^2; M((i-2)*4+3,((i-2)*4)+4) = X(i,1)^3;
    y((i-2)*4+3,1) = Y(i,1);
%Ending slope end condition
M((i-2)*4+4,((i-2)*4)+1) = 0; M((i-2)*4+4,((i-2)*4)+2) = X(i,1)^0; M((i-2)*4+4,((i-2)*4)+3) =2*X(i,1)^1; M((i-2)*4+4,((i-2)*4)+4) = 3*X(i,1)^2;
y((i-2)*4+4,1) = dX(2);
% %Improve matrix rcond
% c = 0.0001;
% M = M + c*eye(size(M));
% rcond(M)
%Evaluate coefficients
a = M\y;
for i = 1:n-1
   output.coefs(i,1) = a(4*i,1); output.coefs(i,2) = a(4*i-1,1); output.coefs(i,3) = a(4*i-2,1); output.coefs(i,4) = a(4*i-3,1); 
end
end

Best Answer

After debugging the code, the solution was very simple than expected. Apparantly I made a mistake in evaluating the cubic spline. The updated code for ppeval_cubic_spline is as follows:
function output = ppval_cubic_spline(pp,sampling_length)
x = [];
y = [];
for i = 1:pp.pieces
    temp_x = pp.breaks(i):sampling_length:pp.breaks(i+1)-sampling_length;
    x = [x temp_x];
    temp_y = pp.coefs(i,1)*temp_x.^3+pp.coefs(i,2)*temp_x.^2+pp.coefs(i,3)*temp_x.^1+pp.coefs(i,4)*temp_x.^0;
    y = [y temp_y];
end
output.x = x;
output.y = y;
end
Related Question