Nvec = [100,500,1000,2000,4000,8000,12000,100000];

for nnn = 1:length(Nvec),
  clear ufe ube;
  N = Nvec(nnn)
  stiff_forced;
  ube = ue - vb;
  ebe(nnn) = sqrt(Tmax/N*ube*ube');
  ufe = ue - vf;
  efe(nnn) = sqrt(Tmax/N*ufe*ufe');
  ute = ue - vt;
  ete(nnn) = sqrt(Tmax/N*ute*ute');
end

hvec = Tmax./Nvec;
loglog(hvec,ebe,'*-');hold on;
loglog(hvec,efe,'*--');
loglog(hvec,ete,'*-.');
grid on;