be(1)=20;ne(1)=20;bf(1)=10;nf(1)=10;

N=20;Qf=0.5063;
for j=1:20;
    r(1)=0;
for t=1:300;
    for i=1:20;
        p(1)=0;
        p(i+1)=p(i)+0.2
     sbe(t)=power(Qf,p(i))*power(ne(t),r(j)*p(i))/(power(20,power(Qf,p(i))*power(ne(t),p(i)*r(j)))+19*power(20,power(Qf,p(i))*power(nf(t),p(i)*r(j))));
     sbf(t)=power(Qf,p(i))*power(nf(t),r(j)*p(i))/(power(20,power(Qf,p(i))*power(ne(t),p(i)*r(j)))+19*power(20,power(Qf,p(i))*power(nf(t),p(i)*r(j))));
    end
     dbe(t)=5*power((1+0.05),t-1)*sbe(t)-0.001*be(t)
     dbf(t)=5*power((1+0.05),t-1)*sbf(t)-0.001*bf(t)
     dne(t)=0.1001*be(t)/1/(1-power(0.1,t))-0.001*ne(t);
     dnf(t)=0.1001*bf(t)/1.2/(1-power(0.1,t))-0.001*nf(t);

     be(t+1)=dbe(t)+ be(t);
     bf(t+1)=dbf(t)+ bf(t);
     ne(t+1)=dne(t)+ne(t);
     nf(t+1)=dnf(t)+nf(t);
end
      r(j+1)=r(j)+0.2;
      be300(j)=be(301);
      bf300(j)=bf(301);
      ne300(j)=ne(301);
      nf300(j)=nf(301);
    x(j)=be300(j)/(be300(j)+19*bf300(j));
    y(j)=bf300(j)/(be300(j)+19*bf300(j));
    A(j)=ne300(j)/(ne300(j)+19*nf300(j));
    B(j)=nf300(j)/(ne300(j)+19*nf300(j));
    HHI(j)=x(j)*x(j)+19*y(j)*y(j);
    HHi(j)=A(j)*A(j)+19*B(j)*B(j);
end
r(21)=[]
plot(r,HHI,'-.',r,HHi, '-*')

我需要把r变成一个常数，而使得be(1)和ne(1)变成一个变量，我应该如何对此程序进行编程，，可以使得我得出be(1)和HHI以及ne(1)和HHi的图形。。。。。。目的就是想分析p的变化对be（1）和ne(1)的变化是如何对HHI产生的影响的。。。
  谢谢谢谢！！！！！！！！！！！！！！！！！！！！！！！