The program

To compute this solution we shall first write a program for freefem to construct a mesh

border(1,0,1,10) begin x:=t; y:=0 end;
border(1,0,1,10) begin x:=1; y:=t end;
border(1,1,0,10) begin x:=t; y:=1 end;
border(1,1,0,10) begin x:=0; y:=t end;
buildmesh(2000);
savemesh("square.msh");

Then we shall write a program which reads the mesh and implement the scheme above. Below the outline of the program is given:

#include <math.h>
#include <stdlib.h>
#include <iostream.h>
#include <fstream.h>
#include <assert.h>
#include "aarray.h"
#include "grid.h"

const int itermax = 100;
void main()
{
	float dt = 0.75/itermax, err = 0, pi = 4*atan(1.0);
	Grid g("square.msh");
	
	A<float> dt2sigma(g.nv), u2(g.nv), u1(g.nv), u0(g.nv);
	
	for(int i=0; i<g.nv;i++)	// initialization
	{	  dt2sigma[i] = 0;
		u0[i] = 0;
		u1[i] = dt*pi*sin(g.v[i].x*pi)*sin(g.v[i].y*pi)*sqrt(2);
	}
	for(int k=0; k<g.nt;k++)
	{  g.t[k].area = ???;
  	 for(int iloc = 0; iloc<3; iloc++)
			    dt2sigma[g(k,iloc)] += g.t[k].area/dt/dt/3;

	for(int iter=0; iter < itermax; iter++)	// time loop
	{
		for(int i=0; i<g.nv;i++)
				u2[i] = 2*u1[i] - u0[i];
		
		for(int k=0; k<g.nt;k++)
			for(int iloc = 0; iloc<3; iloc++)
			 if(!g.t[k].v[iloc]->where)
			 {	int i  = g(k,iloc);
				  int ip = g(k,(iloc+1)%3);
				  int ipp= g(k,(iloc+2)%3);
				  for(int jloc = 0; jloc < 3; jloc++)
				  {
					   int j  = g(k,jloc);
					   float aijk = ???;
					   u2[i] -= aijk * u1[j]/dt2sigma[i];
				  }
			}
		for(int i=0; i<g.nv;i++)
		{  u0[i] = u1[i];
     u1[i] = u2[i];
		}
	}	 
		for(int i=0; i<g.nv;i++)
			 err += pow(u1[i] - sin(itermax*dt*pi*sqrt(2))
                  *sin(g.v[i].x*pi)*sin(g.v[i].y*pi),2);
	 cout << sqrt(err)/g.nv << endl;
}