<< Graphics`PlotField`
<< Graphics`Graphics`

n = 
      30;     
Vinit = 
      0.5;
Vmap = 
        Table[Vinit, {i, 1, n}, {j, 1, 
            n}];

  
Do[Vmap[[1, j]] = 0, {j, 1, n}];
Do[Vmap[[n, j]] = 0, {j, 1, n}];
Do[Vmap[[n/2, j]] = 1, {j, 1, n/2}];
Do[Vmap[[i, 1]] = 1, {i, 1, n}];
Do[
        Vmap[[i, n]] = 0, {i, 1, 
          n}];

  
Niteration = 100;
Do[
     g1 = 
          ListPlot3D[Vmap, PlotRange -> {0, 1}, 
            DisplayFunction -> 
              Identity] ; 
     g2 = 
              ListContourPlot[Vmap, PlotRange -> {0, 1}, 
                DisplayFunction -> 
                  Identity] ;
     Show[
            GraphicsArray[{g1, g2}], 
            DisplayFunction -> $DisplayFunction]; 
Do[
Vmap[[i, 
                j]] = (Vmap[[i - 1, j]] + Vmap[[i + 1, j]] + 
                  Vmap[[i, j - 1]] + Vmap[[i, j + 1]])/
              4, 
{i, 2, n/2 - 1}, {j, 2, 
            n/2}]; 

Do[
Vmap[[i, 
                j]] = (Vmap[[i - 1, j]] + Vmap[[i + 1, j]] + 
                  Vmap[[i, j - 1]] + Vmap[[i, j + 1]])/
              4, 
{i, n/2 + 1, n - 1}, {j, 2, 
            n/2}]; 

Do[
Vmap[[i, 
                j]] = (Vmap[[i - 1, j]] + Vmap[[i + 1, j]] + 
                  Vmap[[i, j - 1]] + Vmap[[i, j + 1]])/
              4, 
{i, 2, n - 1}, {j, n/2 + 1, 
            n - 1}], 
{k, 1, 
          Niteration}];

 (* 
      ListPlot3D[Vmap]; 
ListContourPlot[
        Vmap]; 

Derive FieldMap through  - 
        grad (Vmap[n*n])  
Emap = 
            Table[{Vmap[[i - 1, j]] - 
                  Vmap[[i + 1, j]], -(Vmap[[i, j - 1]] - 
                      Vmap[[i, j + 1]])}, {i, 2, n - 1}, {j, 2, 
                n - 1}]; 
ListPlotVectorField[
        Emap] *)