% Repeating Example 11.1 with Runge function

clear all
close all
clc

f = @(x) 1./(1+25.*x.^2);
a = -1;
b = 1;

% "exact" integral
I = integral(f,a,b);
xplot = linspace(a,b);
yactual = f(xplot);

% First approximation h = 1
h = 1;
x = [a:h:b];
y = f(x);
coeffs = polyfit(x,y,length(x)-1);      % Find the interpolating polynomial

y1plot = polyval(coeffs,xplot);
q = polyint(coeffs);            % Integrate the interpolating polynomial
I1 = diff(polyval(q,[a,b]));    % Evaluate the integrated polynomial at the limits [a,b]

disp('--------------------------------------------------------------------------')
disp(['h = ',num2str(h)])
disp(['Approximate Integral:      ',num2str(I1)])
disp(['Polynomial Coefficients:   ',num2str(coeffs(end:-1:1))])
disp(' ')


figure
hold on
plot(xplot,yactual,'LineWidth',2,'Color',[0, 0.4470, 0.7410],'HandleVisibility','off')
area(xplot,y1plot,'EdgeColor',[0.8500, 0.3250, 0.0980],'FaceAlpha',0.2)
plot(xplot,yactual,'LineWidth',2,'Color',[0, 0.4470, 0.7410])
plot(x,y,'ok','MarkerFaceColor','k','HandleVisibility','off')
xlim([a,b])
legend('Interpolating Polynomial','y_{actual}')
xlabel('x')
ylabel('y')
title(['h = ',num2str(h)])


% Second approximation h = 0.5
h = 0.5;
x = [a:h:b];
y = f(x);
coeffs = polyfit(x,y,length(x)-1);

y1plot = polyval(coeffs,xplot);
q = polyint(coeffs);
I1 = diff(polyval(q,[a,b]));

disp('--------------------------------------------------------------------------')
disp(['h = ',num2str(h)])
disp(['Approximate Integral:      ',num2str(I1)])
disp(['Polynomial Coefficients:   ',num2str(coeffs(end:-1:1))])
disp(' ')

figure
hold on
plot(xplot,yactual,'LineWidth',2,'Color',[0, 0.4470, 0.7410],'HandleVisibility','off')
area(xplot,y1plot,'EdgeColor',[0.8500, 0.3250, 0.0980],'FaceAlpha',0.2)
plot(xplot,yactual,'LineWidth',2,'Color',[0, 0.4470, 0.7410])
plot(x,y,'ok','MarkerFaceColor','k','HandleVisibility','off')
xlim([a,b])
legend('Interpolating Polynomial','y_{actual}')
xlabel('x')
ylabel('y')
title(['h = ',num2str(h)])


% Third approximation h = 0.2
h = 0.2;
x = [a:h:b];
y = f(x);
coeffs = polyfit(x,y,length(x)-1);

y1plot = polyval(coeffs,xplot);
q = polyint(coeffs);
I1 = diff(polyval(q,[a,b]));

disp('--------------------------------------------------------------------------')
disp(['h = ',num2str(h)])
disp(['Approximate Integral:      ',num2str(I1)])
disp(['Polynomial Coefficients:   ',num2str(coeffs(end:-1:1))])
disp(' ')

figure
hold on
plot(xplot,yactual,'LineWidth',2,'Color',[0, 0.4470, 0.7410],'HandleVisibility','off')
area(xplot,y1plot,'EdgeColor',[0.8500, 0.3250, 0.0980],'FaceAlpha',0.2)
plot(xplot,yactual,'LineWidth',2,'Color',[0, 0.4470, 0.7410])
plot(x,y,'ok','MarkerFaceColor','k','HandleVisibility','off')
xlim([a,b])
legend('Interpolating Polynomial','y_{actual}')
xlabel('x')
ylabel('y')
title(['h = ',num2str(h)])


% Fourth approximation h = 0.1
h = 0.1;
x = [a:h:b];
y = f(x);
coeffs = polyfit(x,y,length(x)-1);

y1plot = polyval(coeffs,xplot);
q = polyint(coeffs);
I1 = diff(polyval(q,[a,b]));

disp('--------------------------------------------------------------------------')
disp(['h = ',num2str(h)])
disp(['Approximate Integral:      ',num2str(I1)])
disp(['Polynomial Coefficients:   ',num2str(coeffs(end:-1:1))])
disp(' ')

figure
hold on
plot(xplot,yactual,'LineWidth',2,'Color',[0, 0.4470, 0.7410],'HandleVisibility','off')
area(xplot,y1plot,'EdgeColor',[0.8500, 0.3250, 0.0980],'FaceAlpha',0.2)
plot(xplot,yactual,'LineWidth',2,'Color',[0, 0.4470, 0.7410])
plot(x,y,'ok','MarkerFaceColor','k','HandleVisibility','off')
xlim([a,b])
legend('Interpolating Polynomial','y_{actual}')
xlabel('x')
ylabel('y')
title(['h = ',num2str(h)])