About generalized Stokes problem

Hi, when I try to solve the generallized Stokes problem, I find some Strange phenomena: The pressure error in the cornor is so large. And I don’t know where I was wrong. Who can help me point out the mistakes? Thanks a lot.

The function is

 -2*mu*div(0.5*(grad(u)+grad(u)^T))+grad(ksi)=f
-div(u)+1/lambda*ksi=fp

and the code is

int n = 64;

mesh Th =square(n,n);
//plot(Th,wait=1);

real E = 1000.0, nu = 0.499;
real lambda = E*nu/((1+nu)*(1-2*nu)), mu= E/(2*(1+nu)); //


real sqrt2=sqrt(2.);
macro epsilon(u1,u2) [dx(u1),dy(u2),(dy(u1)+dx(u2))/sqrt2] // EOM
macro div(u,v) ( dx(u)+dy(v) ) // EOM
macro grad(w) [dx(w),dy(w)]  // EOM


// Test case 6
func uxtrue = sin(x);
func uytrue = sin(y);
func ksitrue=cos(x+y);


func fux=2*mu*sin(x)-sin(x+y);
func fuy=2*mu*sin(y)-sin(x+y);
func fp=-(cos(x)+cos(y)+1/lambda*cos(x+y));


fespace Vh(Th,[P2,P2,P1]);
 
    // Define FE variables
Vh [ux,uy,ksi],[vx,vy,phi],[uxerr,uyerr,ksierr];
 


varf GStokes([ux,uy,ksi],[vx,vy,phi]) 
        =int2d(Th)(2*mu*(epsilon(ux,uy)'*epsilon(vx,vy))-ksi*div(vx,vy))
		 +int2d(Th)(fux*vx+fuy*vy)  //rhs 		 
		 +int2d(Th)(-div(ux,uy)*phi-1/lambda*ksi*phi)
		 +int2d(Th)(fp*phi) //rhs
		 +on(1,2,3,4,ux=uxtrue,uy=uytrue);
	
matrix MatB1=GStokes(Vh,Vh,solver=UMFPACK);
real[int] rhs1=GStokes(0,Vh);
ux[]=MatB1^-1*rhs1;

[uxerr,uyerr,ksierr]=[ux,uy,ksi]-[uxtrue,uytrue,ksitrue];
plot(ksi,value=1);
plot(ksierr,value=1);

the Strange places is
image