% Cubic Hermite interpolant for Runge function

f = inline('1 ./ (1 + 25*x.^2)');
npts = 11;
x = [-1:0.2:1]; % 11 equally spaced interpolation pts
y = f(x);% points at which to evaluate interpolant

xx = [-1:.01:1]; % for plotting
fxx = f(xx);

% piecewise cubic Hermite interpolant evaluation
fprime = inline('-50*x ./ (1 + 25*x.^2).^2');
yp = fprime(x);
i = 2;
for ii=1:length(xx)
    while xx(ii) > x(i) & i < npts  % determine which subinterval we are in
        i = i+1;
    end
    h = x(i) - x(i-1);
    alpha = (3/h^2)*(yp(i-1)+yp(i)) + (6/h^3)*(y(i-1)-y(i));
    yy(ii) = - (yp(i-1)/h) * ( (xx(ii)-x(i))^2/2 - h^2/2 ) + ...
        (yp(i)/h) * (xx(ii)-x(i-1))^2/2 + ...
        alpha*(xx(ii)-x(i-1))^2 * ( (xx(ii)-x(i-1))/3 - h/2) + y(i-1);
end

figure(); % creates a new figure window
% plotting the Hermite interpolant 
subplot(2,1,1);
plot(x,y,'o',xx,fxx,'-', xx,yy,'--r');
xlabel('x');
ylabel('1/(1 + 25x^2)');
title('Runge function and its piecewise cubic Hermite interpolant');
legend('Values at nodes','Runge function','Hermite interpolant');

% plotting the errors
subplot(2,1,2);
plot(xx,fxx-yy);
xlabel('x');
ylabel('error');
title('Error of Hermite interpolation');

