function part1gsim1(Llim,numavg)

Lvector = 5:5:Llim;

asympt = [ 1 2 5 14 ];

exactasympt = [];
momsexact = [];

for L=Lvector
  N = L
  c = L/N;
  exact(1) = 1;
  exact(2) = c + 1; 
  exact(3) = c^2 + 3*c + 1 + 1/(N^2);
  exact(4) = c^3 + 6*c^2 + 6*c + 1 + 5*(1 + c)/(N^2);
  
  moms = zeros(1,4);
  for (l=1:numavg) 
    G = (1/sqrt(N))* randn(N,L);
    moms(1) = moms(1) + (1/L)*trace((G'*G));
    moms(2) = moms(2) + (1/L)*trace((G'*G)^2);
    moms(3) = moms(3) + (1/L)*trace((G'*G)^3);
    moms(4) = moms(4) + (1/L)*trace((G'*G)^4);
  end
  moms = moms/numavg;
  
  exactasympt = [ exactasympt (exact - asympt) * (exact - asympt)' ];
  momsexact = [ momsexact (moms - exact) * (moms - exact)' ];
end

plot(Lvector,exactasympt,'k+',Lvector,momsexact,'ko');
legend_handle = legend('MSE between exact and asymptotic moments ','MSE between estimated and exact moments','Location','North');
set(legend_handle,'Box','off');
axis([5 Llim 0 1]);