> restart; with(linalg): with(plots):
> N := 3;
> U := matrix(N+1, N+1); 
> j := 1: for i from 2 to N+1 do U[i,j] := 2.0: od:
> i := 1: for j from 1 to N+1 do U[i,j] := 0.0: od:
> i := N+1: for j from 1 to N+1 do U[i,j] := 1.0: od:
> j := N+1: for i from 2 to N do U[i,j] := 0.0: od:
> print(U);
> for i from 2 to N do
>    for j from 2 to N do
>        eq[i,j] :=U[i,j] =(U[i,j-1] + U[i-1,j] + U[i,j+1] +
> U[i+1,j])/4;
>    od:
> od:
> eq := convert(eq, set);
> solve(eq);
> assign("):  V := submatrix(U, 2..N, 2..N):
> print(V);
> print(U);
> matrixplot(V, axes = framed, style = patchcontour);
> restart; with(linalg): with(plots):
> N := 4;
> U := matrix(N+1, N+1); 
> j := 1: for i from 2 to N+1 do U[i,j] := 2.0: od:
> i := 1: for j from 1 to N+1 do U[i,j] := 0.0: od:
> i := N+1: for j from 1 to N+1 do U[i,j] := 1.0: od:
> j := N+1: for i from 2 to N   do U[i,j] := 0.0: od:
> print(U);
> for i from 2 to N do
>    for j from 2 to N do
>        eq[i,j] :=U[i,j] =(U[i,j-1] + U[i-1,j] + U[i,j+1] +
> U[i+1,j])/4;
>    od:
> od:
> eq := convert(eq, set);
> solve(eq);
> assign("):  V := submatrix(U, 2..N, 2..N):
> print(V);  print(U);
> matrixplot(V, axes = framed, style = patchcontour);
>