Newtonian interpolation polynomial that interpolates f twice (value and 1st derivative)

Construct a Newtonian interpolation polynomial that interpolates f at all points x = −2, −1,0,1,2 twice (value and 1st derivative). What is the value of the interpolation polynomial at x = 1.5?

This is my code but I do not get the correct answer:

format long

function d = deljeneDif_P3_2xInterp(X,Y,DY)
M = zeros(4);
M(:,1) = [Y(1);                      Y(1);  Y(2);                        Y(2)];
M(:,2) = [0;                         DY(1); (Y(2)-Y(1))/(X(2)-X(1));     DY(2)];
M(:,3) = [0;                         0;     (M(3,2)-M(2,2))/(X(2)-X(1)); (M(4,2)-M(3,2))/(X(2)-X(1))];
M(4,4) = (M(4,3)-M(3,3))/(X(2)-X(1));
d = diag(M);
end

a = 1./2.80;

Function definitions in a script must appear at the end of the file.
Move all statements after the "deljeneDif_P3_2xInterp" function definition to before the first local function definition.

Error in connector.internal.fevalMatlab

Error in connector.internal.fevalJSON

f = @(x) (8.*a.^3)./(x.^2+4.*a.^2);
df = @(x) (-16.*a.^3.*x)./(x.^2+4.*a.^2).^2;
interpX = [-2 -1 0 1 2];
interpY = f(interpX);
interpDY = df(interpX);
d = deljeneDif_P3_2xInterp(interpX,interpY,interpDY);
vrednost_2 = polyval(d, 1.5);

Answer 1

[Davide Masiello]

If the problem is the error appearing in your question, i.e.

Function definitions in a script must appear at the end of the file.

then just change the script in the following way and it works.

format long

a = 1./2.80;

f = @(x) (8.*a.^3)./(x.^2+4.*a.^2);
df = @(x) (-16.*a.^3.*x)./(x.^2+4.*a.^2).^2;
interpX = [-2 -1 0 1 2];
interpY = f(interpX);
interpDY = df(interpX);
d = deljeneDif_P3_2xInterp(interpX,interpY,interpDY);
vrednost_2 = polyval(d, 1.5)

vrednost_2 = 
   0.637432906768885

function d = deljeneDif_P3_2xInterp(X,Y,DY)
M = zeros(4);
M(:,1) = [Y(1);                      Y(1);  Y(2);                        Y(2)];
M(:,2) = [0;                         DY(1); (Y(2)-Y(1))/(X(2)-X(1));     DY(2)];
M(:,3) = [0;                         0;     (M(3,2)-M(2,2))/(X(2)-X(1)); (M(4,2)-M(3,2))/(X(2)-X(1))];
M(4,4) = (M(4,3)-M(3,3))/(X(2)-X(1));
d = diag(M);
end

[Nikodin Sedlarevic]

But this value is not the correct answer!