mmOld=Table[1,{x,1,17},{y,1,17}];  (* Set up a table of assumed starting values. Shouldn't matter much *)
(* I originally had 100 grid points, but that was very slow, and didn't improve the look of the final plot *)

(* Set up boundary conditions on the array. Here I'll set it to 2 volts on 2 ends, 0 on the other two sides *)

For[i=1, i ≤17, i++, mmOld[[i,1]]=2];
For[i=1, i≤17, i++, mmOld[[i,17]]=2];
For[j=1,j≤17,j++,mmOld[[1,j]]=0];
For[j=1,j≤17,j++,mmOld[[17,j]]=0];

mmNew=mmOld; 
mdev=1; dev=1;  
(*mdev is the "maximum (squared) change on any single grid point this iteration", while dev will be the sum of the squares of ALL deviations added up *)


For[n=0,dev>.001(* I keep going until the total sum of squared deviations has gotten quite small. Or you could just fix the number of iterations, e.g. n<100 *),n++,  
       mdev=0; dev=0;

(* First sweep through the grid, avoiding the outer row and column, and find the "average" of the old nearest neighbors *)
       For[i=2,i≤16,i++,
             For[j=2,j≤16,j++,  
                   mmNew[[i,j]]= (mmOld[[i-1,j]]+mmOld[[i+1,j]]+mmOld[[i,j-1]]+mmOld[[i,j+1]])/4. ;
                   dev=dev+(mmNew[[i,j]]-mmOld[[i,j]])^2; mdev=Max[mdev,(mmNew[[i,j]]-mmOld[[i,j]])^2]
                   ]
          ];(* For debugging, can Print[{mdev,dev}];*)
(* Now sweep through and update the old grid, and iterate *)  
 For[i=2,i≤16,i++,
        For[j=2,j≤16,j++,
               mmOld[[i,j]]= mmNew[[i,j]] 
               ]
        ]
];
(* Print[{mdev,dev}] *)
ListPlot3D[mmOld, PlotRange -> {0,2}]

(* This is an analytic solution for V=2 on two edges, 0 on the other two *)

v1[x_,y_]:=(8/Pi)*Sum[N[Cosh[(2n+1) Pi *x]] N[Sin[(2n+1) Pi y]]/((2n+1)* N[Cosh[(2n+1) Pi/2]]),{n,0,100}]

Plot3D[v1[xx,yy],{xx,-.5,.5},{yy,0,1}]

(* Now we can set up other boundary conditions .  Here I'll do something like in the homework, where one wall is 2 partway, and then 0 the rest of the way: *)

mmOld=Table[1,{x,1,17},{y,1,17}];
For[i=1, i ≤9, i++, mmOld[[i,1]]=2];
For[i=10, i ≤17, i++, mmOld[[i,1]]=0];
For[i=1, i≤17, i++, mmOld[[i,17]]=0];
For[j=1,j≤17,j++,mmOld[[1,j]]=0];
For[j=1,j≤17,j++,mmOld[[17,j]]=0];
(*mmOld[[1,1]]=1; mmOld[[17,17]]=1.;*)
mmNew=mmOld; 
mdev=1; dev=1;
For[n=0,dev>.001(* n<100 *),n++,  
       mdev=0; dev=0;
       For[i=2,i≤16,i++,
             For[j=2,j≤16,j++,  
                   mmNew[[i,j]]= (mmOld[[i-1,j]]+mmOld[[i+1,j]]+mmOld[[i,j-1]]+mmOld[[i,j+1]])/4. ;
                   dev=dev+(mmNew[[i,j]]-mmOld[[i,j]])^2; mdev=Max[mdev,(mmNew[[i,j]]-mmOld[[i,j]])^2]
                   ]
          ];(* Print[{mdev,dev}];*)  
 For[i=2,i≤16,i++,
        For[j=2,j≤16,j++,
               mmOld[[i,j]]= mmNew[[i,j]] 
               ]
        ]
];

ListPlot3D[mmOld, PlotRange -> {0,2}]

(* One more.  Here I'll do make TWO walls have T=2, but not opposite ones. *)

mmOld=Table[1,{x,1,17},{y,1,17}];
For[i=1, i ≤17, i++, mmOld[[i,1]]=2];
For[i=1, i≤17, i++, mmOld[[i,17]]=0];
For[j=1,j≤17,j++,mmOld[[1,j]]=2];
For[j=1,j≤17,j++,mmOld[[17,j]]=0];
(*mmOld[[1,1]]=1; mmOld[[17,17]]=1.;*)
mmNew=mmOld; 
mdev=1; dev=1;
For[n=0,dev>.001(* n<100 *),n++,  
       mdev=0; dev=0;
       For[i=2,i≤16,i++,
             For[j=2,j≤16,j++,  
                   mmNew[[i,j]]= (mmOld[[i-1,j]]+mmOld[[i+1,j]]+mmOld[[i,j-1]]+mmOld[[i,j+1]])/4. ;
                   dev=dev+(mmNew[[i,j]]-mmOld[[i,j]])^2; mdev=Max[mdev,(mmNew[[i,j]]-mmOld[[i,j]])^2]
                   ]
          ];(* Print[{mdev,dev}];*)  
 For[i=2,i≤16,i++,
        For[j=2,j≤16,j++,
               mmOld[[i,j]]= mmNew[[i,j]] 
               ]
        ]
];

ListPlot3D[mmOld, PlotRange -> {0,2}]

(* This is fun. Let's make some other functional forms*)

mmOld=Table[1,{x,1,17},{y,1,17}];
For[i=1, i ≤17, i++, mmOld[[i,1]]=i/17-1/17];
For[i=1, i≤17, i++, mmOld[[i,17]]=1-1/17 - 2*i/17];
For[j=1,j≤17,j++,mmOld[[1,j]]=-j/17+1/17];
For[j=1,j≤17,j++,mmOld[[17,j]]=-1+1/17];
(*mmOld[[1,1]]=1; mmOld[[17,17]]=1.;*)
mmNew=mmOld; 
mdev=1; dev=1;
For[n=0,dev>.001(* n<100 *),n++,  
       mdev=0; dev=0;
       For[i=2,i≤16,i++,
             For[j=2,j≤16,j++,  
                   mmNew[[i,j]]= (mmOld[[i-1,j]]+mmOld[[i+1,j]]+mmOld[[i,j-1]]+mmOld[[i,j+1]])/4. ;
                   dev=dev+(mmNew[[i,j]]-mmOld[[i,j]])^2; mdev=Max[mdev,(mmNew[[i,j]]-mmOld[[i,j]])^2]
                   ]
          ];(* Print[{mdev,dev}];*)  
 For[i=2,i≤16,i++,
        For[j=2,j≤16,j++,
               mmOld[[i,j]]= mmNew[[i,j]] 
               ]
        ]
];

ListPlot3D[mmOld, PlotRange -> {-1,1}]