% Compute Riemann sums for
% \int_a^b f(x) dx

clear all
hold off

% Define function f(x) 
f = @(x) x.^2
% f = @(x) sqrt(abs(1-x.^2))
% f = @(x) 1./(1+x.^2);
% f = @(x) sin(pi*x);
%f = @(x) tanh(3*x);

% Define interval [a,b]
a=0;
b=1;



% Plot function

nplot = 100;
dxplot = (b-a)/nplot;

xplot= a:dxplot:b;
yplot = f(xplot);

xmin = a;
xmax = b;

ymin = min(yplot);
ymin = min(0,ymin);
ymax = max(yplot);
ymax = max(0,ymax);

plot(xplot,yplot)
axis([xmin,xmax,ymin,ymax]);
xlabel('x');
ylabel('y');


% Loop to compute Riemann sums
% for different numbers of intervals n
while 1 > 0
    % Enter number of intervals n in Riemann sum
    n = input('Enter number of intervals in Riemann sum (0 to end): ');
    if n==0
        break
    end
    n 

    % dx = interval width
    % xl, xr = vectors of left, right x-values
    % yl, yr = vectors of corresponding y-values

    dx = (b-a)/n;
    xl = a:dx:(b-dx);
    yl = f(xl);
    xr = (a+dx):dx:b;
    yr = f(xr);

    

    %Compute Riemann sums

    left_riemann_sum = dx*sum(yl)
    right_riemann_sum = dx*sum(yr)
    average_riemann_sum = 0.5*(right_riemann_sum + left_riemann_sum)

    
    %%% Plot the intervals
    figure(2)
    

subplot(2,1,1)
bar(xr-dx/2,yr,1,'g')
axis([xmin,xmax,ymin,ymax]);

hold on
h = plot(xplot,yplot)
set(h,'LineWidth',4)
set(gca,'LineWidth',2,'FontSize',16)
h = xlabel('x');
set(h,'FontSize',16,'FontWeight','bold')
h = ylabel('y');
set(h,'FontSize',16,'FontWeight','bold')
h = title(strcat('Right Endpoint Sum = ', num2str(right_riemann_sum)))
set(h,'FontSize',16,'FontWeight','bold')
hold off


subplot(2,1,2)
bar(xl+dx/2,yl,1,'g')
axis([xmin,xmax,ymin,ymax]);

hold on
h = plot(xplot,yplot)
set(h,'LineWidth',4)
set(gca,'LineWidth',2,'FontSize',16)
h = xlabel(strcat('Average Riemann Sum = ',num2str(average_riemann_sum)))
set(h,'FontSize',16,'FontWeight','bold')
h = ylabel('y');
set(h,'FontSize',16,'FontWeight','bold')
h = title(strcat('Left Endpoint Sum = ', num2str(left_riemann_sum)))
set(h,'FontSize',16,'FontWeight','bold')

hold off


    
end