Don't post code as an image.

The original post contains code as an image.

Thanks Jim. Yes, I figured it that the integration limits were to be changed in order to avoid this issue. The code works now!

Yes, sure,.....thanks for your response. I am sure you'll have a better way to code this same problem, but I am quite a beginner with MATHEMATICA. I have added print commands at every output, which is my way to understand as I code,.....

Clear["Global`*"]

(* Parameters *)
Em = 1;
v = 0.2;
rm = 2;

(* a - Deriving the 3 Shape Functions *)
(* Create Node Locations *)
posN1 = {2, 0};
posN2 = {5, 0};
posN3 = {2, 2};

(* Displacement Approximation of the Element*)
DisplacementU = a0 + a1*X1 + a2*X2;
DisplacementV = b0 + b1*X1 + b2*X2;

(* Writing the Approximation functions in terms of Nodal DOFs for the \
Element *)
 DisplacementUN1 = a0 + a1*posN1[[1]] + a2*posN1[[2]] ;
 DisplacementUN2 = a0 + a1*posN2[[1]] + a2*posN2[[2]] ;
 DisplacementUN3 = a0 + a1*posN3[[1]] + a2*posN3[[2]] ;
 DisplacementVN1 = b0 + b1*posN1[[1]] + b2*posN1[[2]] ;
 DisplacementVN2 = b0 + b1*posN2[[1]] + b2*posN2[[2]] ;
 DisplacementVN3 = b0 + b1*posN3[[1]] + b2*posN3[[2]] ;

(* Solving for the Coefficients in Terms of the Nodal Displacements & \
Nodal Positions, for that Element *)
sol1 = Solve[{DisplacementUN1 == u1, DisplacementUN2 == u2, 
    DisplacementUN3 == u3}, {a0, a1, a2}];
{a0, a1, a2} = {a0, a1, a2} /. sol1[[1]];
Print["  "];
Print["The {a0,a1,a2} Coefficients Are"];
Print[{a0, a1, a2}];
sol2 = Solve[{DisplacementVN1 == v1, DisplacementVN2 == v2, 
    DisplacementVN3 == v3}, {b0, b1, b2}];
{b0, b1, b2} = {b0, b1, b2} /. sol2[[1]];
Print["  "];
Print["The {b0,b1,b2} Coefficients Are"];
Print[{b0, b1, b2}];

(* The solved Coefficients fit into the Displacement Approximation \
Function for the Element *)
Print["  "];
Print["The X1-direction Displacement Approximation Function for the \
Element is"];
Print[DisplacementU];
Print["  "];
Print["The X2-direction Displacement Approximation Function for the \
Element is"];
Print[DisplacementV];

(* However, the 'Displacement Approximation Function' has to be \
re-arranged to obtain the Shape Functions for that Element *)

(* Solving the Shape Functions for the Element *)
N1 = Coefficient[DisplacementU, u1];
N2 = Coefficient[DisplacementU, u2];
N3 = Coefficient[DisplacementU, u3];
Ni = {N1, N2, N3} ;
Print["  "];
Print["The Shape Functions Are "];
Print[ Ni];
Print["  "];

(* Create Strain-Displacement Matrix [B] of the Element *)
B = {{D[N1, X1], 0, D[N2, X1], 0, D[N3, X1], 0}, {N1/X1, 0, N2/X1, 0, 
    N3/X1, 0}, {0, D[N1, X2], 0, D[N2, X2], 0, D[N3, X2]}, {D[N1, X2],
     D[N1, X1], D[N2, X2], D[N2, X1], D[N3, X2], D[N3, X1]}};
Print["The [B] Matrix is "];
Print[B // MatrixForm];
Print["  "];

(* Create the Constitutive Matrix [Cc] of the Element *)
Cc = ((Em*(1 - v))/((1 + v)*(1 - 2*v)))*{{1, v/(1 - v), v/(1 - v), 
     0}, {v/(1 - v), 1, v/(1 - v), 0}, {v/(1 - v), v/(1 - v), 1, 
     0}, {0, 0, 0, (1 - 2*v)/(2*(1 - v))}};
Print["The Constitutive Matrix [Cc] of the Element is "];
Print[Cc // MatrixForm];
Print["  "];

(* Integration Function For Stiffness Matrix*)
IntFunc = (Transpose[B] . Cc . B)*(rm + X1);

(* Stiffness Matrix - Exact Integration *)
sol3 = NIntegrate[IntFunc, {X1, 0, 3}, {X2, 0, 2}, AccuracyGoal -> 1, 
   WorkingPrecision -> 10, PrecisionGoal -> 2];
Print["The Stiffness Matrix [Ke] of the Element is "];
Print[sol3 // MatrixForm];