clear all;clc;format compact;format long;

% matrix of a, b, c coefficients:

% rows: Fe, S, Fe1, Mo

% columns: A1; B1; A2; B2; A3; B3; A4; B4; C

fM = ...

[11.9185 4.87394 7.04848 0.34023 3.34326 15.9330 2.27228 79.0339 1.40818;...

 7.18742 1.43280 5.88671 0.02865 5.15858 22.1101 1.64403 55.4651 -3.87732;...

 11.9185 4.87394 7.04848 0.34023 3.34326 15.9330 2.27228 79.0339 1.40818;...

 19.3885 0.97877 11.8308 10.0885 3.75919 31.9738 1.46772 117.932 5.55047];

%%% store radius data:

% distance from: origin

% columns: Fe, S, Fe, Mo

R_el = [2.0 3.3 3.5 3.5];

RHO_t = zeros(4,400);

 for numel = 1:4 

 EL = numel;

 RHO = zeros(1,400);

 dmax = zeros(1,400);

     for iter = 1:400

        dmax(iter) = iter/100; % in angstroms

% numerical integration

     int_fun = @(s) 4*pi*(s.^2).* ...

    (fM(EL,1).*exp(-fM(EL,2).*(s.^2)*0.25) + ...

     fM(EL,3).*exp(-fM(EL,4).*(s.^2)*0.25) + ...

     fM(EL,5).*exp(-fM(EL,6).*(s.^2)*0.25) + ...

     fM(EL,7).*exp(-fM(EL,8).*(s.^2)*0.25) + fM(EL,9)).* ...

     sin(2*pi*s*R_el(EL))./(2*pi*s*R_el(EL));

     RHO(iter) = quad(int_fun,0,1/dmax(iter));

        clc;display(iter);display(numel);

     end

 RHO_t(numel,:) = RHO;

 end

RHO_t(1,:)= 6*RHO_t(1,:);

RHO_t(2,:)= 9*RHO_t(2,:);

 figure;

 axis([0.5 3.5 -10 10]); hold on;

 plot(dmax,RHO_t(1,:),...

      dmax,RHO_t(2,:),...

      dmax,RHO_t(3,:),...

      dmax,RHO_t(4,:),...

      dmax,sum(RHO_t,1));

  title('Electron Density Profile');

  legend('Fe','S','Fe1','Mo','Sum');

  xlabel('d_m_a_x'); ylabel('Rho(r)');

  set(gca,'XDir','reverse');

##############

Any suggestions will be appreciated. Thanks! 

Niu