zz=1:82;
subplot(2,2,1);
  plot(zz,zz.^2);hold on;
  plot(zz,1/4*zz.^2)
  plot(zz,1/9*zz.^2)
  plot(zz,1/16*zz.^2) 
set(gca,'ylim',[0 85]) 
set(gca,'fontsiz',14)
line([2 2],[0 100],'color',[0.7 0.7 0.7])
line([10 10  ],[0 100],'color',[0.7 0.7 0.7])
line([18 18],[0 100],'color',[0.7 0.7 0.7])
line([36 36 ],[0 100],'color',[0.7 0.7 0.7])
title('Z^2 / n^2')
 

subplot(2,2,2);
  plot(zz,zz.^2);hold on;
  plot(zz,1/4*zz.^2)
  plot(zz,1/9*zz.^2)
  plot(zz,1/16*zz.^2)
set(gca,'ylim',[0 85]) 
set(gca,'fontsiz',14)
  yy1=1:2;plot(yy1,yy1.^2,'y*')
  ztot=yy1.^2;
  ytot= yy1;
  yy2=3:10;plot(yy2,1/4*yy2.^2,'y*')
ytot=[ytot yy2 ];
ztot=[ztot 1/4*yy2.^2];

  yy3=11:18;plot(yy3,1/9*yy3.^2,'y*')
ytot=[ytot yy3 ];
ztot=[ztot 1/9*yy3.^2];
  yy4=19:36;plot(yy,1/16*yy.^2,'y*')
ytot=[ytot yy4 ];
ztot=[ztot 1/16*yy4.^2];
plot(ytot,ztot,'k-','linewidth',2)
set(gca,'fontsiz',14)
line([2 2],[0 100],'color',[0.7 0.7 0.7])
line([10 10  ],[0 100],'color',[0.7 0.7 0.7])
line([18 18],[0 100],'color',[0.7 0.7 0.7])
line([36 36 ],[0 100],'color',[0.7 0.7 0.7])
title('Z^2 / n^2  with shell filling')

subplot(2,2,3);
  plot(zz,zz.^2,':');hold on;
  plot(zz,1/4*(zz-2).^2,':')
  plot(zz,1/9*(zz-10).^2,':')
  plot(zz,1/16*(zz-18).^2,':')

  yy1=1:2;plot(yy1,yy1.^2,'y*')
  ztot=yy1.^2;
  ytot= yy1;
  yy2=3:10;plot(yy2,1/4*(yy2-2).^2,'y*')
ytot=[ytot yy2 ];
ztot=[ztot 1/4*(yy2-2).^2];

  yy3=11:18;plot(yy3,1/9*(yy3-10).^2,'y*')
ytot=[ytot yy3 ];
ztot=[ztot 1/9*(yy3-10).^2];
  yy4=19:36;plot(yy,1/16*(yy4-18).^2,'y*')
ytot=[ytot yy4 ];
ztot=[ztot 1/16*(yy4-18).^2];
plot(ytot,ztot,'k-','linewidth',2)
set(gca,'ylim',[0 20]) 
set(gca,'fontsiz',14)
line([2 2],[0 100],'color',[0.7 0.7 0.7])
line([10 10  ],[0 100],'color',[0.7 0.7 0.7])
line([18 18],[0 100],'color',[0.7 0.7 0.7])
line([36 36 ],[0 100],'color',[0.7 0.7 0.7])

title('(Z-n_s)^2 / n^2  with shell filling')
set(gca,'fontsiz',14)

subplot(2,2,4);

  yy1=1:2;plot(yy1,yy1.^2,'y*')
  ztot=yy1.^2;
  ytot= yy1;
  yy2=3:4; 
ytot=[ytot yy2 ];
ztot=[ztot 1/4*(yy2-2).^2];
  yy22=5:6; 
ytot=[ytot yy22 ];
ztot=[ztot 1/4*(yy22-4).^2];
  yy23=7:10; 
ytot=[ytot yy23 ];
ztot=[ztot 1/4*(yy23-6).^2];


  yy3=11:12;plot(yy3,1/9*(yy3-10).^2,'y*')
ytot=[ytot yy3 ];
ztot=[ztot 1/9*(yy3-10).^2];

  yy3=13:14;plot(yy3,1/9*(yy3-12).^2,'y*')
ytot=[ytot yy3 ];
ztot=[ztot 1/9*(yy3-12).^2];

  yy3=15:18;plot(yy3,1/9*(yy3-14).^2,'y*')
ytot=[ytot yy3 ];
ztot=[ztot 1/9*(yy3-14).^2];


  yy4=19:20;
ytot=[ytot yy4 ];
ztot=[ztot 1/16*(yy4-18).^2];

  yy4=21:24;
ytot=[ytot yy4 ];
ztot=[ztot 1/16*(yy4-20).^2];

  yy4=25:28;
ytot=[ytot yy4 ];
ztot=[ztot 1/16*(yy4-24).^2];

  yy4=29:36;
ytot=[ytot yy4 ];
ztot=[ztot 1/16*(yy4-30).^2];

ytot=[ytot 37 ];
ztot=[ztot 0.5];

plot(ytot,ztot,'k-','linewidth',2)

set(gca,'ylim',[0 5]) 
set(gca,'fontsiz',14)
line([2 2],[0 100],'color',[0.7 0.7 0.7])
line([10 10  ],[0 100],'color',[0.7 0.7 0.7])
line([18 18],[0 100],'color',[0.7 0.7 0.7])
line([36 36 ],[0 100],'color',[0.7 0.7 0.7])

title('(Z- n_v )^2 / n^2  - variable n_v and with shell filling')
set(gca,'fontsiz',14)
     %%%
 set(gcf,'color','white')