Simpson's Rule Program, modifying it to print intervals when new values are added

I'm using the following program to approximate an integral of a function using Simpson's Rule. I want to modify it so I can print the intervals for when new values are added to simpson_result and also print values of f(x) over the entire interval [a,b] so I can sketch the curves on paper. I'm not sure where to print these values and also since I input the function as an argument in the command window I'm not sure how to edit the program to print values of f(x). Thanks for any help!

function simpson_result=recur_simpson(Fun,a,b,TOL,level,level_max);
%%function simpson_result=recur_simpson(Fun,a,b,TOL,level,level_max);
% INPUT
% Fun - function
% a,b - lower and upper limits
% TOL - desired accuracy of the integral
% level - level of recursion (set=0)
% level_max - maximum depth of recursion (set=4)
% OUTPUT
% simpson_result - value of the integral
level=level+1;
h=b-a;
c=(b+a)/2;
one_simpson=h*(feval(Fun,a)+4*feval(Fun,c)+feval(Fun,b))/6;
d=(a+c)/2;
e=(c+b)/2;
two_simpson=h*(feval(Fun,a)+4*feval(Fun,d)+2*feval(Fun,c)+4*feval(Fun,e)+feval(Fun,b))/12;
if level >= level_max 
    simpson_result=two_simpson;
    message='Level Limit Exceeded'
    elseif abs(two_simpson-one_simpson)<15*TOL;
        simpson_result=two_simpson+(two_simpson-one_simpson)/15;
    else
        left_simpson=recur_simpson(Fun,a,c,TOL/2,level,level_max);
        right_simpson=recur_simpson(Fun,c,b,TOL/2,level,level_max);
        simpson_result=left_simpson+right_simpson;
end;