% Example Rectangle Method
clearvars
close all


% Set spacing h
h = 0.2;

f = @(x) x.^2;
a = 0;
b = 1;
x1 = [a:h:b];
x2 = repelem(x1(2:end-1),2);
xh = [a x2 b];

% Exact Integral
I = integral(f,a,b);
xplot = linspace(a,b);
yactual = f(xplot);

% Numerical Integrals
I1=0;
I2=0;
I3=0;
for i=1:length(x1)-1
    I1 = I1+f(a+(i-1)*h)*h;
    I2 = I2+f(a+(i-1/2)*h)*h;
    I3 = I3+f(a+(i)*h)*h;
end
e1 = I-I1;
e2 = I-I2;
e3 = I-I3;

% Rectangle Method 1
for i=1:length(xh)
    if mod(i,2)==0
        y1(i)=f(xh(i-1));
    else
        y1(i)=f(xh(i));
    end
end

% Rectangle Method 2
for i=1:length(xh)
    if mod(i,2)==0
        y2(i)=(f(xh(i))+f(xh(i-1)))/2;
    else
        y2(i)=(f(xh(i+1))+f(xh(i)))/2;
    end
end

% Rectangle Method 3
for i=1:length(xh)
    if mod(i,2)==0
        y3(i)=f(xh(i));
    else
        y3(i)=f(xh(i+1));
    end
end

% Plotting
figure(1)
plot(xplot,yactual,'LineWidth',2)
hold on
area(xh,y1,'EdgeColor',[0.8500, 0.3250, 0.0980],'FaceAlpha',0.2)
xlabel('x')
ylabel('y')
title(['I1 = ',num2str(I1)])
hold off

figure(2)
plot(xplot,yactual,'LineWidth',2)
hold on
area(xh,y2,'EdgeColor',[0.8500, 0.3250, 0.0980],'FaceAlpha',0.2)
xlabel('x')
ylabel('y')
title(['I2 = ',num2str(I2)])
hold off

figure(3)
plot(xplot,yactual,'LineWidth',2)
hold on
area(xh,y3,'EdgeColor',[0.8500, 0.3250, 0.0980],'FaceAlpha',0.2)
xlabel('x')
ylabel('y')
title(['I3 = ',num2str(I3)])
hold off