```
load "Element_P3"
load "Element_P4"
int k = 1*2*2;
mesh Th = square(10*k, 10*k);
//plot(Th);
func Pk = P4;
fespace Vhc ( Th, Pk );
Vhc uc,ucX,ucY,fC,fCh,ucEX,err;
//--exact solution, partial derivatives and rhs----
ucEX = exp(x^2/4+y^2/4);
ucX = exp(x^2/4+y^2/4)*x/2;
ucY = exp(x^2/4+y^2/4)*y/2;
fC = exp(x^2/4+y^2/4)*(-x^2/4-y^2/4);
//------ get the matrice and right hand side ------
varf ahc(u,v) = int2d(Th)( dx(u)*dx(v)+dy(u)*dy(v) +u*v );
varf rhsC(unused,v) = int2d(Th)(fC*v)
+int1d(Th,1)( (N.x*ucX+N.y*ucY)*v )
+int1d(Th,2)( (N.x*ucX+N.y*ucY)*v )
+int1d(Th,3)( (N.x*ucX+N.y*ucY)*v )
+int1d(Th,4)( (N.x*ucX+N.y*ucY)*v );
//----------------solve the linear system----------
matrix Ahc;
Ahc = ahc(Vhc,Vhc);
set(Ahc, solver=CG, eps=1e-20);
fCh[] = rhsC(0,Vhc);
uc[] = Ahc^-1*fCh[];
//--------------- get the error -------------------
err[] = uc[]-ucEX[];
cout.scientific<<"Infinity error:"
<<max(abs(err[].min), abs(err[].max))<<endl;
real l2er;
l2er = ( int2d(Th)( (ucEX-uc)^2) );
cout.scientific<<"L2 error: "<<sqrt(l2er)<<endl;
```

**Question: in this example, I verified the L2 convergence rate for P3 element as k varies which is correct.**

**However, the convergence rate for P4 element is incorrect. Am I wrong somewhere?**

Thanks a lot for your help.

P4Element.edp (1.1 KB)`Preformatted text`