(*technique constants*)
alpha=1;
damp=1;
x0=-1;
x1=1;
y0=-1;
y1=1;
t0=0;
t1=5;
dx=.1;
dy=.1;
dt=.01;

(*derived constants*)
(*Courant numbers*)
Cx=alpha*dt/dx;
Cy=alpha*dt/dy;

(*number of points*)
Dx=IntegerPart[(x1-x0)/dx]+1;
Dy=IntegerPart[(y1-y0)/dy]+1;
timeLength=(t1-t0)/dt+1;

(*derived functions*)

partitionToDomain[{a_,b_}]:=N[{x0+(a-1)*dx,y0+(b-1)*dy}];
domainToPartition[{x_,y_}]:=IntegerPart[{(x-x0)/dx+1,(y-y0)/dy+1}];
indexToPartition[n_]:=IntegerPart[{Mod[n-1,Dx]+1,(n-(Mod[n-1,Dx]+1))/Dx+1}];
partitionToIndex[{a_,b_}]:=IntegerPart[(b-1)*Dx+a];

indexOfTime[t_]:=IntegerPart[(t-t0)/dt -1];

(*state vector*)
u0=Table[0,{i,1,Dx*Dy}];
b=Table[0,{i,1,Dx*Dy}];
(*initial conditions*)
(*f[x_,y_]:=E^(-10x^2)*E^(-10y^2);*)
f[x_,y_]:=0;
(*t=t0 time boundary conditions*)
For[x=-1,x≤1,x+=dx,
For[y=-1,y≤1,y+=dy,
u0[[partitionToIndex[domainToPartition[{x,y}]]]]=Chop[N[f[x,y]],.001];
];
];
For[x=-.0,x≤0,x+=dx,
For[y=-0,y≤0,y+=dy,
u0[[partitionToIndex[domainToPartition[{x,y}]]]]=1;
];
];

(*spatial boundary conditions*)

dirichlet=IdentityMatrix[Dx*Dy];
(*along the bottom of the partition*)
For[i=1,i≤Dx,i++,
dirichlet[[partitionToIndex[{i,1}]]]=ConstantArray[0,Dx*Dy];
];
(*along the top of the partition*)
For[i=1,i≤Dx,i++,
dirichlet[[partitionToIndex[{i,Dy}]]]=ConstantArray[0,Dx*Dy];
];
(*along the left of the partition*)
For[i=1,i≤Dy,i++,
dirichlet[[partitionToIndex[{1,i}]]]=ConstantArray[0,Dx*Dy];
];
(*along the right of the partition*)
For[i=1,i≤Dy,i++,
dirichlet[[partitionToIndex[{Dx,i}]]]=ConstantArray[0,Dx*Dy];
];

(*setup matrix*)
(*Interior points*)
A=SparseArray[{Band[{1,1}]->(2-2Cx^2-2Cy^2),Band[{1,2}]->Cx^2,Band[{2,1}]->(Cx^2),Band[{1,1+Dx}]->(Cy^2),Band[{1+Dx,1}]->(Cy^2)},{Dx*Dy,Dx*Dy}];
(*edge points*)
(*edge indices
bottom edge 1 to Dx (exclusive)
left edge k*Dx+1 for k from 0 to Dy-1 (exclusive)
right edge k*Dx for k from 0 to Dy (exclusive)
top edge Dx*(Dy-1)+1 to Dx*Dy (exclusive)*)
For[a=2,a<Dx,a++,
A[[partitionToIndex[{a,1}]]]=ConstantArray[0,Dx*Dy];
A[[partitionToIndex[{a,1}]]][[partitionToIndex[{a,1}]]]=2-2Cx^2-Cy^2;
A[[partitionToIndex[{a,1}]]][[partitionToIndex[{a,1}]-1]]=Cx^2;
A[[partitionToIndex[{a,1}]]][[partitionToIndex[{a,1}]+1]]=Cx^2;
A[[partitionToIndex[{a,1}]]][[partitionToIndex[{a,1}]+Dx]]=Cy^2;
];
For[k=2,k<Dy,k++,
A[[partitionToIndex[{1,k}]]]=ConstantArray[0,Dx*Dy];
A[[partitionToIndex[{1,k}]]][[partitionToIndex[{1,k}]]]=2-Cx^2-2Cy^2;
A[[partitionToIndex[{1,k}]]][[partitionToIndex[{1,k}]-Dx]]=Cy^2;
A[[partitionToIndex[{1,k}]]][[partitionToIndex[{1,k}]+Dx]]=Cy^2;
A[[partitionToIndex[{1,k}]]][[partitionToIndex[{1,k}]+1]]=Cx^2;
];
For[k=2,k<Dy,k++,
A[[partitionToIndex[{1,k}]]]=ConstantArray[0,Dx*Dy];
A[[partitionToIndex[{1,k}]]][[partitionToIndex[{1,k}]]]=2-Cx^2-2Cy^2;
A[[partitionToIndex[{1,k}]]][[partitionToIndex[{1,k}]-Dx]]=Cy^2;
A[[partitionToIndex[{1,k}]]][[partitionToIndex[{1,k}]+Dx]]=Cy^2;
A[[partitionToIndex[{1,k}]]][[partitionToIndex[{1,k}]-1]]=Cx^2;
];
For[a=2,a<Dx,a++,
A[[partitionToIndex[{a,Dy}]]]=ConstantArray[0,Dx*Dy];
A[[partitionToIndex[{a,Dy}]]][[partitionToIndex[{a,Dy}]]]=2-2Cx^2-Cy^2;
A[[partitionToIndex[{a,Dy}]]][[partitionToIndex[{a,Dy}]-1]]=Cx^2;
A[[partitionToIndex[{a,Dy}]]][[partitionToIndex[{a,Dy}]+1]]=Cx^2;
A[[partitionToIndex[{a,Dy}]]][[partitionToIndex[{a,Dy}]-Dx]]=Cy^2;
];
(*corner points*)
(*corner indices are 1(bottom left), Dx(bottom right), Dx*(Dy-1)+1(upper left), and Dx*Dy(upper right)*)
A[[partitionToIndex[{1,1}]]]=ConstantArray[0,Dx*Dy];
A[[partitionToIndex[{1,1}]]][[partitionToIndex[{1,1}]]]=2-Cx^2-Cy^2;
A[[partitionToIndex[{1,1}]]][[partitionToIndex[{1,1}]+1]]=Cx^2;
A[[partitionToIndex[{1,1}]]][[partitionToIndex[{1,1}]+Dx]]=Cy^2;

A[[partitionToIndex[{Dx,1}]]]=ConstantArray[0,Dx*Dy];
A[[partitionToIndex[{Dx,1}]]][[partitionToIndex[{Dx,1}]]]=2-Cx^2-Cy^2;
A[[partitionToIndex[{Dx,1}]]][[partitionToIndex[{Dx,1}]-1]]=Cx^2;
A[[partitionToIndex[{Dx,1}]]][[partitionToIndex[{Dx,1}]+Dx]]=Cy^2;

A[[partitionToIndex[{1,Dy}]]]=ConstantArray[0,Dx*Dy];
A[[partitionToIndex[{1,Dy}]]][[partitionToIndex[{1,Dy}]]]=2-Cx^2-Cy^2;
A[[partitionToIndex[{1,Dy}]]][[partitionToIndex[{1,Dy}]+1]]=Cx^2;
A[[partitionToIndex[{1,Dy}]]][[partitionToIndex[{1,Dy}]-Dx]]=Cy^2;

A[[partitionToIndex[{Dx,Dy}]]]=ConstantArray[0,Dx*Dy];
A[[partitionToIndex[{Dx,Dy}]]][[partitionToIndex[{Dx,Dy}]]]=2-Cx^2-Cy^2;
A[[partitionToIndex[{Dx,Dy}]]][[partitionToIndex[{Dx,Dy}]-1]]=Cx^2;
A[[partitionToIndex[{Dx,Dy}]]][[partitionToIndex[{Dx,Dy}]-Dx]]=Cy^2;

(*initial kick*)
u1=Table[0,{i,1,Dx*Dy}];

g[x_,y_]:=0;
For[x=-1,x≤1,x+=dx,
For[y=-1,y≤1,y+=dy,
u1[[partitionToIndex[domainToPartition[{x,y}]]]]=Chop[N[g[x,y]],.001];
];
];
u1[[partitionToIndex[domainToPartition[{0,0}]]]]=0;
u1=u0+dt*u1+.5*A.u0+.5*b;

(*solutionbuilding loop*)
answer={u0};
u2=ConstantArray[0,Dx*Dy];
For[t=t0,t≤t1,t+=dt,
u2=Chop[(1/(1+2*damp*dt))(A.u1+2*damp*dt*u1-u0+b),.001];
u0=u1;
u1=u2;
answer=Join[answer,{u2}];
];

timeslices=Table[Table[{partitionToDomain[indexToPartition[i]][[1]],partitionToDomain[indexToPartition[i]][[2]],answer[[t]][[i]]},{i,1,Dx*Dy}],{t,1,Length[answer]}];

skipfactor=1;
rasterized=Table[Image[ListPlot3D[timeslices[[skipfactor*i]],PlotRange->{-1,2},Mesh->None,InterpolationOrder->1]],{i,1,IntegerPart[Length[timeslices]/skipfactor]-1}];
Export["2dwaves.avi",rasterized,"FrameRate"->60]

(*skipfactor=1;
rasterized=Table[Image[ListPlot3D[timeslices[[skipfactor*i]],PlotRange->{-1,2},ColorFunction->ColorData[{"NeonColors","Reverse"}],Mesh->None,InterpolationOrder->0]],{i,1,IntegerPart[Length[timeslices]/skipfactor]-1}];
Export["neonhd.gif",rasterized,"FrameRate"->60]*)