Just remove the undesirable equations and variables:

Clear@"`*";
diffFormula = (T[i + 1, k] - 2 T[i, k] + T[i - 1, k])/x^2 + 
              (T[i, k + 1] - 2 T[i, k] + T[i, k - 1])/y^2 == 0;

NP = 5;(*number of grid points*)
x = y = 0.05;
Table[T[0, k] = 325., {k, 0, NP}];
Table[T[i, NP] = 325., {i, 0, NP}];
Table[T[i, 0] = 325., {i, 0, NP}];
Table[T[NP, k] = 325., {k, 0, NP}];

inner = {T[2, 2], T[2, 3], T[3, 2], T[3, 3]};

eq = Flatten[Table[diffFormula, {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedeq = Select[eq, Count[#, Alternatives @@ inner, Infinity] < 2 &];

var = Flatten[Table[T[i, k], {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedvar = Complement[var, inner];

(# = 400) & /@ inner;

solution = Flatten@Solve[selectedeq, selectedvar];
s = Table[T[i, k], {k, 0, NP}, {i, 0, NP}] /. solution;

ListPlot@s

BTW, though the region you're dealing with is luckily simple, handling irregular region with FDM can be really frustrating. This post is an example. If FDM isn't necessary for you, have a look at the FEM capabilities of NDSolve new added in v10.

Just remove the undesirable equations and variables:

Clear@"`*";
diffFormula = (T[i + 1, k] - 2 T[i, k] + T[i - 1, k])/x^2 + 
              (T[i, k + 1] - 2 T[i, k] + T[i, k - 1])/y^2 == 0;

NP = 5;(*number of grid points*)
x = y = 0.05;
Table[T[0, k] = 325., {k, 0, NP}];
Table[T[i, NP] = 325., {i, 0, NP}];
Table[T[i, 0] = 325., {i, 0, NP}];
Table[T[NP, k] = 325., {k, 0, NP}];

inner = {T[2, 2], T[2, 3], T[3, 2], T[3, 3]};

eq = Flatten[Table[diffFormula, {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedeq = Select[eq, Count[#, Alternatives @@ inner, Infinity] < 2 &];

var = Flatten[Table[T[i, k], {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedvar = Complement[var, inner];

(# = 400) & /@ inner;

solution = Flatten@Solve[selectedeq, selectedvar];
s = Table[T[i, k], {k, 0, NP}, {i, 0, NP}] /. solution;

ListPlot@s

BTW, though the region you're dealing with is luckily simple, handling irregular region with FDM can be really frustrating. This post is an example. If FDM isn't necessary for you, have a look at the FEM capabilities of NDSolve new added in v10.

Just remove the undesirable equations and variables:

Clear@"`*";
diffFormula = (T[i + 1, k] - 2 T[i, k] + T[i - 1, k])/x^2 + 
              (T[i, k + 1] - 2 T[i, k] + T[i, k - 1])/y^2 == 0;

NP = 5;(*number of grid points*)
x = y = 0.05;
Table[T[0, k] = 325., {k, 0, NP}];
Table[T[i, NP] = 325., {i, 0, NP}];
Table[T[i, 0] = 325., {i, 0, NP}];
Table[T[NP, k] = 325., {k, 0, NP}];

inner = {T[2, 2], T[2, 3], T[3, 2], T[3, 3]};

eq = Flatten[Table[diffFormula, {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedeq = Select[eq, Count[#, Alternatives @@ inner, Infinity] < 2 &];

var = Flatten[Table[T[i, k], {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedvar = Complement[var, inner];

(# = 400) & /@ inner;

solution = Flatten@Solve[selectedeq, selectedvar];
s = Table[T[i, k], {k, 0, NP}, {i, 0, NP}] /. solution;

ListPlot@s

BTW, handling irregular region with FDM can be really frustrating. This post is an example.

Just remove the undesirable equations and variables:

Clear@"`*";
diffFormula = (T[i + 1, k] - 2 T[i, k] + T[i - 1, k])/x^2 + 
              (T[i, k + 1] - 2 T[i, k] + T[i, k - 1])/y^2 == 0;

NP = 5;(*number of grid points*)
x = y = 0.05;
Table[T[0, k] = 325., {k, 0, NP}];
Table[T[i, NP] = 325., {i, 0, NP}];
Table[T[i, 0] = 325., {i, 0, NP}];
Table[T[NP, k] = 325., {k, 0, NP}];

inner = {T[2, 2], T[2, 3], T[3, 2], T[3, 3]};

eq = Flatten[Table[diffFormula, {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedeq = Select[eq, Count[#, Alternatives @@ inner, Infinity] < 2 &];

var = Flatten[Table[T[i, k], {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedvar = Complement[var, inner];

(# = 400) & /@ inner;

solution = Flatten@Solve[selectedeq, selectedvar];
s = Table[T[i, k], {k, 0, NP}, {i, 0, NP}] /. solution;

ListPlot@s

BTW, though the region you're dealing with is luckily simple, handling irregular region with FDM can be really frustrating. This post is an example. If FDM isn't necessary for you, have a look at the FEM capabilities of NDSolve new added in v10.