function [Qout,fcount] = simp_grdemo(f,a,b,tol)
% implements composite SIMPSON'S rule by halving the step size repeatedly;
% shows the steps in approximating the integral:
%   the integral over the subinterval currently being computed is shaded,
%   and the color switches to green when the desired tolerance is reached.

%--------------------------------------------------------------
shg
clf reset
set(gcf,'doublebuffer','on','menubar','none', ...
   'numbertitle','off','name','Simpson graphical demo')
%--------------------------------------------------------------

% initial evaluations of f
c = (a + b)/2;
fa = f(a);
fc = f(c);
fb = f(b);

%--------------------------------------------------------------
h = (b-a)/2;
x = [a c b];
y = [fa fc fb];
maxy = max(y);
miny = min(y);
for k = 1:31
   v = f(a+k*h/32);
   maxy = max(maxy,v);
   miny = min(miny,v);
end
set(gcf,'userdata',0)
hold on
p(1) = fill(a,fa,'k');
p(2) = fill(b,fb,'k');
p(3) = plot(x,y,'.','markersize',16);
hold off
s = (maxy - miny)/20;
axis([a b miny-s maxy+s])
q(1) = uicontrol('string','step', ...
   'units','normal','pos',[.65 .02 .08 .06], ...
   'callback','set(gcf,''userdata'',1)');
q(2) = uicontrol('string','auto', ...
   'units','normal','pos',[.75 .02 .08 .06], ...
   'callback','set(gcf,''userdata'',2)');
q(3) = uicontrol('string','quit', ...
   'units','normal','pos',[.85 .02 .08 .06], ...
   'callback','set(gcf,''userdata'',3)');
%--------------------------------------------------------------

% recursive call to subfunction that divides interval in two 
% and applies rule to each half
[Q,k] = simp_grstep(f,a,b, tol, fa,fc,fb);
fcount = 3 + k;

%--------------------------------------------------------------
title(sprintf('Q = %8.4f,   fcount = %4.0f',Q,fcount))
delete(p(1:2));
delete(q(1:2));
set(q(3),'string','close','callback','close(gcf)')
if nargout > 0, Qout = Q; end
%--------------------------------------------------------------

end

% subfunction

function [Q,k] = simp_grstep(f,a,b,tol,fa,fc,fb)

h=(b-a)/2; 
c=(a+b)/2;
% divide into two subintervals and evaluate f at new nodes
d=(a+c)/2;
e=(c+b)/2;
fd = f(d);
fe = f(e);
% apply quadrature formula
Qprev = (1/3)*h*(fa + 4*fc + fb);
Qnew = (1/3)*h/2*(fa + 4*fd + 2*fc + 4*fe + fb);

%--------------------------------------------------------------
u1 = a:h/32:c;
v1 = polyinterp([a d c],[fa fd fc],u1);
u1 = [a u1 c];
v1 = [0 v1 0];
u2 = c:h/32:b;
v2 = polyinterp([c e b],[fc fe fb],u2);
u2 = [c u2 b];
v2 = [0 v2 0];
if abs(Qnew - Qprev) <= tol
   color = [0 2/3 0];
else
   color = [.6 .6 .6];
end
p = flipud(get(gca,'child'));
x = [get(p(3),'xdata') d e];
y = [get(p(3),'ydata') fd fe];
set(p(1),'xdata',u1,'ydata',v1,'facecolor',color)
set(p(2),'xdata',u2,'ydata',v2,'facecolor',color)
set(p(3),'xdata',x,'ydata',y)
set(gca,'xtick',sort(x),'xticklabel',[]);
title(num2str(length(x)))
pause(.25)
while get(gcf,'userdata') == 0
   pause(.25)
end
if get(gcf,'userdata') == 1
   set(gcf,'userdata',0)
end
%--------------------------------------------------------------

% check convergence
if (abs(Qnew - Qprev) <= tol) | (get(gcf,'userdata') == 3)
   Q  = Qnew;
   k = 2;
else  % repeat the subdivision process
   [Ql,kl] = simp_grstep(f,a,c, tol/2, fa,fd,fc);
   [Qr,kr] = simp_grstep(f,c,b, tol/2, fc,fe,fb);
   Q  = Ql + Qr;
   k = 2 + kl + kr;
end

end
