Use plot(x,y,'o',xx,yy) to plot the given points together with the interpolating curve (the vectors x, y contain the coordinates of the given points, the vectors xx, yy contain the coordinates of the points on the interpolating curve).

function y = f1(x)
global nf            % declare nf as global variable
y = x.*cos(x).^3;    % must use ".*","./",".^" since x may be vector
nf = nf + length(x); % update total number of function evaluations

Then you can use quad and find the number of function evaluations as follows:

global nf           % declare nf as global (IMPORTANT!!)
nf = 0              % reset nf
q = quad('f1',0,1,.001)
nf                  % show number of function evaluations

Note that you have to reset nf before each call of quad.

Note that a straight forward implementation is inefficient since it re-evaluates the function at the same points several times. Once you've got the straight forward implementation working, rewrite it as follows: Use the following adquad.m

function q = adquad(func,a,b,tol)
fa = feval(func,a); fb = feval(func,b);
q = adquad1(func,a,b,tol,fa,fb);

and write a function adquad1 which only uses one function evaluation fc=feval(func,c) at the midpoint c. If the estimated error in adquad1 is too large, it calls itself again for the two subintervals, passing on the already computed function values in the last two arguments.