Romberg code isn't running correctly.

I am trying to write code for Romberg integration and am getting a error message when I run the program. Can someone help me with the problem please. Thanks

   function Rnum = MSipesR(f,a,b,m)
  % Integration algorithm based on Romberg extrapolation
  % f - string input for function y = f(x) (e.g. f = 'x.^6')
  % a - lower limit of integration
  % b - upper limit of integration
  % m - maximal number of Romberg iterations
  % Rnum - row-vector of numerical approximations for the integral of f(x) between x = a and x = b:
  % the entry with index k in Rnum corresponds to the approximation of order O(h^(2k))

R = ones(m,m); % the matrix for Romberg approximations 

hmin = (b-a)/2^(m-1); % the minimal step size 

for k = 1 : m

      h = 2^(k-1)*hmin; % the step size for the k-th row of the Romberg matrix

  x = a : h : b;
  f1 = eval(f);
  k1 = length(f1);
    R(k,1) = 0.5*h*(f1(1) + 2*sum(f1(2:k1-1)) + f1(k1)); % trapezoidal rule for the first column of the Romberg matrix
end

for k = 2 : m % compute the k-th column of the Romberg matrix

      for kk = 1 : (m-k+1) % the matrix of Romberg approximations is triangular!

          R(kk,k) = R(kk,k-1)+(R(kk,k-1) - R(kk+1,k-1))/(4^(k-1)-1); % see the Romberg integration algorithm
      end 

end

% define the vector Rnum for numerical approximations
Rnum = R(1,:);

[Are Mjaavatten]

What is the problem? It seems to work very nicely for me:

 Rnum = MSipesR('x.^6',0,1,4)
 Rnum =
    0.1506    0.1430    0.1429    0.1429

Check, using quad:

 quad(@(x) x.^6,0,1)
 ans =
    0.1429