f=zetac(x)

function f=zetac(x)
% f=zetac(x) computes the Riemann zeta function
% for real or complex argument

% Riemann zeta function from dissertation of
%    Professer Robert W. Scharstein. 
% The function handles scalar or matrix arguments.
%           RWS, 13 June 2005
if nargin==0
  [u,v]=meshgrid(linspace(0,1,40),linspace(5,30,80)); 
  x=u+i*v;
end
siz=size(x); x=x(:); M=20; q=1; n=0:M;
x2=x.*(x+1).*(x+2);
x4=x2.*(x+3).*(x+4);
x6=x4.*(x+5).*(x+6);
first=zeros(size(x)); nq=n+q;
for k=1:M+1; first=first+nq(k).^(-x); end
f=first+1./((x-1).*(M+q).^(x-1))-(1/2)./(M+q).^x...
    +x./(12*(M+q).^(x+1))-x2./(720*(M+q).^(x+3))...
    +x4./(30240*(M+q).^(x+5))-x6./(1209600*(M+q).^(x+7)); 
f=reshape(f,siz); 
if nargin>0, return, end
surf(u,v,abs(f)),view([77,16]), xlabel('x -axis')
ylabel('y axis'), title('abs( zeta( x + i*y ) )');
f=[]; shg