MATLAB index exceeds matrix dimensions error.

I have created the following script to specifically implement piecewise quadratic interpolation. The script seems to work for several inputs however for the examples shown below I recieve an error: index exceeds matrix dimensions.

SCRIPT:

 function [p,px,pxx] = quadinterp(x, xj, fj)

 n = length(xj);

 [~,j] = histc(x,xj); 

 j(x>=xj(n)) = n-1;  

 j(x<xj(1)) = 1; 

 p = fj(j).*((x - xj(j+1)).*(x-xj(j+2)))./((xj(j) - xj(j+1)).*(xj(j) - xj(j+2))) ...
    + fj(j+1).*((x - xj(j)).*(x-xj(j+2)))./((xj(j+1) - xj(j)).*(xj(j+1) - xj(j+2)))...
    + fj(j+2).*((x - xj(j)).*(x-xj(j+1)))./((xj(j+2) - xj(j)).*(xj(j+2) - xj(j+1)));

 px = fj(j).*(2*x-xj(j+1)-xj(j+2))./((xj(j) - xj(j+1)).*(xj(j) - xj(j+2))) ...
    + fj(j+1).*(2*x-xj(j)-xj(j+2))./((xj(j+1) - xj(j)).*(xj(j+1) - xj(j+2)))...
    + fj(j+2).*((2*x-xj(j)-xj(j+1)))./((xj(j+2) - xj(j)).*(xj(j+2) - xj(j+1)));

 pxx = fj(j).*2./((xj(j) - xj(j+1)).*(xj(j) - xj(j+2))) ...
    + fj(j+1).*2./((xj(j+1) - xj(j)).*(xj(j+1) - xj(j+2)))...
    + fj(j+2).*(2)./((xj(j+2) - xj(j)).*(xj(j+2) - xj(j+1)));

The following inputs give correct solutions:

 xj = [1 5 9 11 13];

 fj = [xj].^2;

 x = linspace(2,10,13);

The following outputs give the error:

 xj = [1 5 9 11 13];

 fj = [xj].^2;

 x = linspace(1,12,13);
&
 xj = linspace(0, 2,5); 

 fj = xj.^2;

 x = linspace(0, 2, 20);

Any help is greatly appreciated.

Answer 1

[Azzi Abdelmalek]

Check this line

   p = fj(j).*((x - xj(j+1)).*(x-xj(j+2)))./((xj(j) - xj(j+1)).*(xj(j) - xj(j+2))) ...
    + fj(j+1).*((x - xj(j)).*(x-xj(j+2)))./((xj(j+1) - xj(j)).*(xj(j+1) - xj(j+2)))...
    + fj(j+2).*((x - xj(j)).*(x-xj(j+1)))./((xj(j+2) - xj(j)).*(xj(j+2) - xj(j+1)))

you have

j=[1 1 1 1 1 2 2 2 2 3 3 4 4]

the max value of j is 4, and the size of xj is 5, when you type xj(j+2), you are asking for x([3 3 3 3 3 4 4 4 4 5 5 6 6]), x(6) can't be calculated, because the size of xj is [1 6]

[JLK]

Thanks for your input.

Is there anyway of correcting this without changing the line

 p = fj(j).*((x - xj(j+1)).*(x-xj(j+2)))./((xj(j) - xj(j+1)).*(xj(j) - xj(j+2))) ...
    + fj(j+1).*((x - xj(j)).*(x-xj(j+2)))./((xj(j+1) - xj(j)).*(xj(j+1) - xj(j+2)))...
    + fj(j+2).*((x - xj(j)).*(x-xj(j+1)))./((xj(j+2) - xj(j)).*(xj(j+2) - xj(j+1)))

as this is the formula for lagrange interpolation