(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



InternalEq:= (
k=i+1;
A[[k,k-1]]=  1.0;
A[[k,k      ]]=-2.0;
A[[k,k+1]]=   1.0;
)
InternalRHS:=(
k=i+1;
b[[k]]=f[a+i*\[CapitalDelta]x,n*\[CapitalDelta]t]
)
RobinBoundaryEq:=(
k=i+1;
Switch[ i,
0,   mult=  1; x=a,
nx, mult=-1; x=b
,_, Print["Error"]
];
A[[k,k]]=-2.0+(1 + mult \[Alpha] \[CapitalDelta]x);
A[[k,k+mult]]= 1.0
)
RobinRHS:=(
k=i+1;
Switch[ i,
0,   mult=   1; x=a;
nx, mult=-1; x =b,
_, Print["Error"]
];
b[[k]]=f[x,n*\[CapitalDelta]t]- mult  \[CapitalDelta]x g[x,n*\[CapitalDelta]t];
)
PeriodicBoundaryEq:=(
k=i+1;
Switch[ i,
0, x=a,
nx,x=b,
_, Print["Error"]
];
A[[k,k]]=-2.0;
A[[k,2]]=  1.0;
A[[k,nx-1]]=1.0;
)
PeriodicRHS:=(
k=i+1;
Switch[ i,
0,   mult=   1; x=a,
nx, mult=-1; x =b,
_, Print["Error"]
];
b[[k]]=f[x,n*\[CapitalDelta]t];
)
IC:=(
k=i+1;
U0[[k]]=u0[a+i*\[CapitalDelta]x];
)


InternalEq[{i_,j_},{\[Nu]x_,\[Nu]y_},{nx_,ny_},k_]:= {
{k[i,j],k[i,j]}->1-2 \[Nu]x-2\[Nu]y,
{k[i,j],k[i-1,j]}      ->\[Nu]x,
{k[i,j],k[i+1,j]}      ->\[Nu]x,
{k[i,j],k[i,      j-1]}->\[Nu]y,
{k[i,j],k[i,      j+1]}->\[Nu]y
}


PeriodicEq[{i_,j_},{\[Nu]x_,\[Nu]y_},{nx_,ny_},k_]:={
{k[i,j],k[i,j]}->1-2 \[Nu]x-2\[Nu]y,
{k[i,j],k[Mod[i-1,nx],j]}      ->\[Nu]x,
{k[i,j],k[Mod[i+1,nx],j]}      ->\[Nu]x,
{k[i,j],k[i,      Mod[j-1,ny]]}->\[Nu]y,
{k[i,j],k[i,      Mod[j+1,ny]]}->\[Nu]y
}


uxNeumannEq[{i_,j_},{\[Nu]x_,\[Nu]y_},{nx_,ny_},k_]:={
{k[i,j],k[i,j]}->1-\[Nu]x-2\[Nu]y,
Evaluate[
Switch[i,
0,   {k[i,j],k[i+1,j]}->\[Nu]x,
nx, {k[i,j],k[i-1,j]}->\[Nu]x]
],
{k[i,j],k[i,      j-1]}->\[Nu]y,
{k[i,j],k[i,      j+1]}->\[Nu]y
}
uyNeumannEq[{i_,j_},{\[Nu]x_,\[Nu]y_},{nx_,ny_},k_]:=Module[
{},
{
{k[i,j],k[i,j]}->1-2\[Nu]x-\[Nu]y,
 {k[i,j],k[i+1,j]}->\[Nu]x,
{k[i,j],k[i-1,j]}->\[Nu]x,
Switch[j,
0,  {k[i,j],k[i,j+1]}->\[Nu]y,
ny,{k[i,j],k[i,j-1]}->\[Nu]y
]
}
]


CornerEq[{i_,j_},{\[Nu]x_,\[Nu]y_},{nx_,ny_},k_]:={
{k[i,j],k[i,j]}->1-2\[Nu]x-\[Nu]y,
{k[i,j],k[Mod[i+1,nx],j]}->\[Nu]x,
{k[i,j],k[Mod[i-1,nx],j]}->\[Nu]x,
Switch[j,
0,  {k[i,j],k[i,j+1]}->\[Nu]y,
ny,{k[i,j],k[i,j-1]}->\[Nu]y
]
}


A[{\[Nu]x_,\[Nu]y_},{nx_,ny_}]:=Module[
{i,j,k,SysSize=(nx+1)*(ny+1)},
k[i_,j_]:= 1+i+j*(nx+1);
SparseArray[
Flatten[
{
Table[  InternalEq[{i,j},{\[Nu]x,\[Nu]y},{nx,ny},k], 
 {i,1,nx-1,1},   {j,1,ny-1,1}],
Table[  PeriodicEq[{i,j},{\[Nu]x,\[Nu]y},{nx,ny},k], 
 {i,1,nx-1,1},   {j,0,ny,      ny}],
Table[  uyNeumannEq[{i,j},{\[Nu]x,\[Nu]y},{nx,ny},k],
  {i,0,nx,       nx},{j,1,ny-1,1}]
(*,
Table[   CornerEq[{i,j},{\[Nu]x,\[Nu]y},{nx,ny},k],
  {i,0,nx,       nx},{j,0,ny      ,ny}] *)
}
],
{SysSize, SysSize}]
]



