I am trying to solve this integration by simpsons method and plot figure. It's give figure taking only the value of P0= -6 from for loop. How can I solve it?

   alpha = 45;     
  beta = 185;     

gamma_e = 116;  

% Gain values

G_ei = -18.96;    
G_ee = 18.52;
G_sr = -0.26;
G_rs = 16.92;
G_es = 2.55;
G_re = 4.67;
G_se = 0.73;
G_sn = 2.78;

G_esre = G_es*G_sr*G_re;
G_srs = G_sr*G_rs;
G_ese = G_es*G_se;
G_esn = G_es*G_sn;

t_0 = 0.085;    % corticothalamic loop delay in second

r_e = 0.086;    % Excitatory axon range in metre

f = linspace(-40,40,500);   % f = frequency in Hz

w = 2*pi*f;                % angular frequency in radian per second

delt_P = 0.5;

L=zeros(1,500);
Q=repmat(L,1);
P=repmat(L,1);

%%%%%%%%%%%%%integration %%%%%%%%%%%%

a = -80*pi;
b = 80*pi;
n=500;

I=repmat(L,1);
P_initial = repmat(L,1);
P_shift = repmat(L,1);
p = repmat(L,1);

for k = 1:length(w)
    for P0 = [6 -6]

L_initial = @(w1) (1-((1i*w1)/alpha))^(-1)*(1-((1i*w1)/beta))^(-1);                                                                 

Q_initial = @(w1) (1/(r_e^2))*((1-((1i*w1)/gamma_e))^(2) - (1/(1-G_ei*L_initial(w1)))*....
    (L_initial(w1)*G_ee + (exp(1i*w1*t_0)*(L_initial(w1)^2*G_ese +L_initial(w1)^3*G_esre))/(1-L_initial(w1)^2*G_srs)));                  

P_initial =  @(w1) (pi/r_e^4)* (abs((L_initial(w1)^2*G_esn)/((1-L_initial(w1)^2*G_srs)*....
    (1-G_ei*L_initial(w1)))))^2 * abs((atan2((imag(Q_initial(w1))),(real(Q_initial(w1)))))/imag(Q_initial(w1)));                   

G =  150*exp(- (f - P0).^2./(2*(delt_P).^2));   

P2 = @(w1) G(k) + P_initial(w1);

L_shift =  @(w1) (1-((1i*(w(k)-w1))/alpha))^(-1)* (1-((1i*(w(k)-w1))/beta))^(-1);                               

Q_shift  = @(w1)  (1/(r_e^2))*((1-((1i*(w(k)-w1))/gamma_e))^(2) - (1/(1-G_ei*L_shift(w1)))*...
    (L_shift(w1)*G_ee + (exp(1i*(w(k)-w1)*t_0)*(L_shift(w1)^2*G_ese +L_shift(w1)^3*G_esre))/(1-L_shift(w1)^2*G_srs)));     

P_shift =  @(w1)  (pi/r_e^4)* (abs((L_shift(w1)^2*G_esn)/((1-L_shift(w1)^2*G_srs)*(1-G_ei*L_shift(w1)))))^2 *....
    abs((atan2((imag(Q_shift(w1))),(real(Q_shift(w1)))))/imag(Q_shift(w1)));                                

p =  @(w1)  P2(w1)*P_shift(w1);        %  Power spectrum formula for P(w1)*p(w-w1)

I(k) = simprl(p,a,b,n);

      end
  end

figure(1)
plot(f,I,'r--')

figure(2)
plot(f,G,'k')

%%%%% Simprl %%%

    function [s] = simprl(p,a,b,n) % f is the function to be integrated; a = initial value of the interval; b = final value of the interval; n = No. of subintervals 
    % The function implements the Simpson's Rule

    h = (b-a)./n;

    s1 = 0;  %  The variable s1 is initialised to 0.
    s2=0;    %  The variable s2 is initialised to 0.

    % loop for odd values in the range

    for k = 1:n/2;              % The index variable k starts at 1, then increases in steps of 1 until it reaches n/2.
        x = a + h*(2*k-1);
        s1 = s1+feval(p,x);      % Each time through the loop the value of feval(f,x) is added to s1.
    end

    % loop for even values in the range

    for k = 1:(n/2 - 1);
        x = a + h*2*k;
        s2 = s2+feval(p,x);   
    end

    % Final result of integration where odd values are multiplied by 4 and even values are multiplied by 2

    s = h*(feval(p,a)+feval(p,b)+4*s1+2*s2)/3;

[Walter Roberson]

http://www.mathworks.co.uk/matlabcentral/answers/43073-a-guide-to-tags

Answer 1

[Walter Roberson]

You have the structure

for k = something
  for P0 = something
    ...
    I(k) = simprl(p,a,b,n);
  end
end

Each time you iterate through the different P0's, you store only one value, the value for the active k. The effect is as if you had only done the calculation for the last P0.

You should be indexing your I according to k and also according to how many P0 you have done.

I(k, P0idx) = ...

where P0idx would be something you had set at 1 for the first P0 index, at 2 for the next P0 index, and so on.

[Namira]

Thanks. I faced difficulties to set index of P0. Would you please suggest me how can I set it?

[Walter Roberson]

P0vals = [6 -6];
for P0idx = 1 : length(P0vals)
  P0 = P0vals(P0idx);
  ...
  I(k, P0idx) = ...
end

[Namira]

Thanks. It helps me a lot.