clear figures
modelLB1 = createpde();

square = [3 4 0 1 1 0 1 1 0 0]';
g = decsg(square);
geometryFromEdges(modelLB1,g);

subplot(2,2,1);
pdegplot(modelLB1,'EdgeLabels','on');
axis equal

applyBoundaryCondition(modelLB1,'Edge',1,'u',@(region,state)region.x+1);
applyBoundaryCondition(modelLB1,'Edge',2,'u',@(region,state)1+region.y);
applyBoundaryCondition(modelLB1,'Edge',3,'u',@(region,state)region.x);
applyBoundaryCondition(modelLB1,'Edge',4,'u',@(region,state)region.y);

subplot(2,2,2);
generateMesh(modelLB1,'Hmax',0.1,'GeometricOrder','quadratic'); %quadratic
pdemesh(modelLB1);
axis equal

specifyCoefficients(modelLB1,'m',0,'d',0,'c',1,'a',0,'f',0);
                      
u = solvepde(modelLB1); 
u = u.NodalSolution;

subplot(2,2,3);
pdeplot(modelLB1, 'XYData',u)
title('Solution'); xlabel('x'); ylabel('y');
axis equal

subplot(2,2,4);
p = modelLB1.Mesh.Nodes;
exact = p(1,:) + p(2,:);
pdeplot(modelLB1,'XYData',u - exact')
title('Error'); xlabel('x'); ylabel('y')
axis equal