function shells(fstr, a, b)

    func = str2func(sprintf('@(x) %s', fstr));

    % Create a figure and axes
    fig = figure('Visible','off');
    %ax = axes('Units','pixels');
    
    % Create slider
    notches = 50 
    sld = uicontrol('Style', 'slider',...
        'Min',1,'Max',notches+1,'Value',2,...
        'Position', [100 20 300 20],...
        'SliderStep', [1 5]/notches,...
        'Callback', @reslide); 

    % Add a text uicontrol to label the slider.
    txt = uicontrol('Style','text',...
        'Position',[400 20 120 20],...
        'String','N');
    
    plotpow(2);
    campos([-3 -4 2]);

    %surf(peaks)
    % Make figure visble after adding all components
    fig.Visible = 'on';
    % This code uses dot notation to set properties. 
    % Dot notation runs in R2014b and later.
    % For R2014a and earlier: set(f,'Visible','on');

    function reslide(source, callbackdata)
        val = source.Value;
        plotpow(val);
    end

    function plotpow(N)
        N = ceil(N);
        txt.String = sprintf('N = %d', N);
        k = 2;
        x = linspace(a, b, N+1);
        cla;
        volume = zeros(size(x));
        for i = 1:N
            %cap(x, x.^k);
            h = 0.5 * (func(x(i)) + func(x(i+1)));
            shell(x(i), x(i+1), h, h);
            volume(i) = pi * (x(i+1)^2 - x(i)^2) * h;
            %annulus(x(i), x(i+1), f(x(i)), f(x(i+1)))
            %annulus(x(i), x(i+1), x(i)^k, x(i+1)^k)
        end
        title(sprintf('Volume: %.6f', sum(volume)));
        axis equal
    end

    function cap(x0, radius)
        [r, theta] = meshgrid(linspace(0, radius, 5), linspace(0, 2*pi, 41));
        hold on
        surf(x0 * ones(size(r)), r .* cos(theta), r .* sin(theta));
        hold off
    end

    function annulus(x0, x1, h0, h1)
        [h, theta] = meshgrid([h0 h1], linspace(0, 2*pi, 41)); 
        [r, theta] = meshgrid([x0 x1], linspace(0, 2*pi, 41)); 
        hold on
        surf(r .* cos(theta), r .* sin(theta), h);
        hold off
    end

    function shell(x0, x1, h0, h1)
        annulus(x0, x1, 0, 0);
        annulus(x0, x0, 0, h0);
        annulus(x1, x1, 0, h1);
        annulus(x0, x1, h0, h1);
    end

    function cylinder(x0, x1, r0)
        cap(x0, r0);
        annulus(x0, x1, r0, r0);
        cap(x1, r0);
    end
end
