em1=[];em2=[];inc=0.01;em1a1=[];em1m1=[];em2a1=[];em2m1=[];
yav=[];yng=[];beta=[];adn=[];add=[];gamma=[];res=[];
E=15;L=0.020;R=50;T=1/10;
for m_1=1:1:10
   m_2=m_1*10;
mu=0:inc:1;
   em1(1,:)=2.71.^(m_1*(1-mu(1,:))*T);
   em2(1,:)=2.71.^(m_2*(1-mu(1,:))*T);
beta(m_1,:)=(m_2*em1(1,:)-m_1*em2(1,:))/(m_2-m_1);
gamma(1,:)=(em1(1,:)-em2(1,:))/(m_2-m_1);
em1a1(1,:)=(1+em1(1,:));
em1m1(1,:)=(1-em1(1,:));
em2a1(1,:)=(1+em2(1,:));
em2m1(1,:)=(1-em2(1,:));
coef1=E/L*T; coef2=E/L; coef3=E/R;
adn(1,:)=m_2*em1a1(1,:)-m_1*em2a1(1,:);
add(1,:)= m_2*em1m1(1,:)-m_1*em2m1(1,:);
yav(m_1,:)=( adn(1,:).*mu(1,:)*0.5*coef1 )./add(1,:)+...
   +( (em2(1,:)-em1(1,:))*E/L )./add(1,:)+...
   +(E/R)*( m_2-m_1-m_2*em1(1,:)+m_1*em2(1,:) )./add(1,:);
yng(m_1,:)=((m_2*em1(1,:)-m_1*em2(1,:)).*mu(1,:)*coef1)/add(1,:)+...
      +( (em2(1,:)-em1(1,:))*coef2 )./add(1,:)+...
      +(E/R)*( m_2-m_1-m_2*em1(1,:)+m_1*em2(1,:) )./add(1,:);
end
plot(yav(1,:),yng(1,:),yav(2,:),yng(2,:),yav(3,:),yng(3,:),yav(4,:),yng(4,:),...
   yav(5,:),yng(5,:),yav(6,:),yng(6,:),yav(7,:),yng(7,:),yav(8,:),yng(8,:),...
   yav(9,:),yng(9,:),yav(10,:),yng(10,:));
axis([1 2e04 1 1600]);
plot(yng(1,:),yav(1,:),yng(2,:),yav(2,:),yng(3,:),yav(3,:),yng(4,:),yav(4,:),...
   yng(5,:),yav(5,:),yng(6,:),yav(6,:),yng(7,:),yav(7,:),yng(8,:),yav(8,:),...
   yng(9,:),yav(9,:),yng(10,:),yav(10,:));
axis([1 1600 1 16e3]);
title('y_{average} vs y^{-}_{\infty} several root positive')
xlabel('y^{-}_{\infty}')
ylabel(' y_{average}')
