Numerical integration in plane polar space

I accidentally posted this in the newsgroup but it is better posted here

I am writing a program to work out the helmholtz free energy between two uniaxial birefringent materials. This involves integrating a function, log(D0(n,r,phi)) over the k-space r=0:infinity and phi=0:2pi using the limit on r, r=0:1e21. I have tried to weight each segment by it's area A=r x dphi x dr, however this causes my outputs to jump in magnitude past their expacted values, and also causes erratic behaviour of the function. Below is the code, where the functions defining D0 are condensed into integrand(n,r,phi,theta)

~~~~~~ function omega=helmholtz(theta)

D_ncounter=zeros(N+1,1);

for n=0:1:N D_rcounter=zeros(N+1,length(rspace)+1);

for rr=1:length(rspace) r=rspace(rr);

             D0(n+1,rr,pp)=integrand(n,r,phi,theta);

D_phicounter(n+1,rr,pp)=(R/length(rspace))*r*(2*pi/length(phispace))*log(D0(n+1,rr,pp)); %helmholtz free energy for one dr and dphi (multiplied by the weight of each step)

        end %end of phi loop

          D_rcounter(n+1,rr) = D_rcounter(n+1,rr) + sum(D_phicounter(n+1,rr,:)); %helmholtz free energy for one dr, summing over all dphis (multiplied by the weight of each ring of r+-dr)
            %D_rcounter(n+1,rr)
    end %end of r loop

    if n==0
        D_ncounter(n+1,1)=D_ncounter(n+1,1)+0.5*sum(D_rcounter(n+1,:));
    else
        D_ncounter(n+1,1)=D_ncounter(n+1,1)+sum(D_rcounter(n+1,:));
    end %end of if

end %end of n loop

omega=sum(D_ncounter(:,1)); ~~~~~~

The results of this code, for a range of fixed thetas, are put into an array, M, and differentiated using diff(M)/dphi, is this the correct method of differentiating the results?

Thanks in advance