in the example code, the source term is

prode = -0.126*kp*(pow(2*dx(Ux),2)+pow(2*dy(Uy),2)+2*pow(dx(Uy)+dy(Ux),2))/2;
prodk = -prode*kp/epp*0.09/0.126;

but both of them should be negative, so prodk should be

prodk = prode*kp/epp*0.09/0.126;

in the example code, the source term is

prode = -0.126*kp*(pow(2*dx(Ux),2)+pow(2*dy(Uy),2)+2*pow(dx(Uy)+dy(Ux),2))/2;
prodk = -prode*kp/epp*0.09/0.126;

but both of them should be negative, so prodk should be

prodk = prode*kp/epp*0.09/0.126;