natural cubic spline interpolation

I need to write a function that takes a list of x values and an equations as parameters and outputs a 4xn matrix made out of the coefficients of the relative cubic polynomials. (e.g. the kth row of the matrix represents the cubic interpolation between the points and )

I wrote (read: coopted) this function but it is giving an error on line 29 (coefmat(i,:)=[s0,s1,s2,s3];) saying that the matrices have different dimensions

function coefmat = splinecalc(Xlist, f)
    Ylist=zeros(1,length(Xlist));
    coefmat=zeros(length(Xlist)-1,4);
    for i=1:length(Xlist)
        Ylist(1,i)=f(Xlist(1,i));
    end
    for i=1:length(Xlist)
        n = length(Xlist);
        h = Xlist(2:n) - Xlist(1:n-1);
        d = (Ylist(2:n) - Ylist(1:n-1))./h;
        lower = h(1:end-1);
        main  = 2*(h(1:end-1) + h(2:end));
        upper = h(2:end);

        T = spdiags([lower main upper], [-1 0 1], n-2, n-2);
        rhs = 6*(d(2:end)-d(1:end-1));

        m = T\rhs;

        % Use natural boundary conditions where second derivative
        % is zero at the endpoints

        m = [ 0; m; 0];

        s0 = Ylist(1,i);
        s1 = d - h.*(2*m(1:end-1) + m(2:end))/6;
        s2 = m/2;
        s3 =(m(2:end)-m(1:end-1))./(6*h);
        coefmat(i,:)=[s0,s1,s2,s3];
    end
    coefmat;
end

Answer 1

[John D'Errico]

Learn to use the debugger. This is a basic problem. You know what line the problem is on, but even if you did not, then you could set the debugger to stop when it hits an error.

dbstop if error

But simplest is to just set a breakpoint at that line. Then you would ask two questions of yourself.

What is the size of coefmat? Hint:

coefmat=zeros(length(Xlist)-1,4);

coefmat has 4 columns.

What are these variables: s0,s1,s2,s3?

Pay particular attention to s3, which would seem to be a vector of numbers itself. If there is a bug in this code however, you might want to find better code that lacks bugs.

So think about the size of coefmat. It seems to be n array with 4 columns. What does the offending line of code do?

coefmat(i,:)=[s0,s1,s2,s3];

It tries to stuff some numbers into a row of coefmat. So it needs a 1x4 vector of numbers. Then look at what you are trying to put in there.

s0 seems to be a scalar.

s1 seems to be a scalar.

s2 seems to be a scalar.

s3 seems to be a vector, of length 1 less than the length of m.

Can you stuff more than 4 elements into a vector of length 4? (NO.)

Sorry, but I won't debug this buggy code, to tell you what it should do, because that would force me to verify ALL of the essentially uncommented code.

Find better code if you cannot write it yourself. Anyway, why do you think you need to use a natural cubic spline, anyway? You would usually be better off using that which spline itself produces. That is, spline also gives you a cubic spline, but with a better choice of end conditions than the natural ones, which are often an issue themselves.