(* ::Package:: *)

(*initialisation.m*)

sd=10^6; tol=10^(-5); nt=1; 

(*Provide N, the "resolution" of the system we wish to simulate*)
NStep=6;
iini=0;ilast=NStep;
jini=0;jlast=NStep;
(*End of provide N, the "resolution" of the system we wish to simulate*)

(*The coordinates of the boundaries*)
x[iini]=-1.0;x[ilast]=1.0;
y[jini]=-1.0;y[jlast]=1.0;
Deltax=(x[ilast]-x[iini])/NStep;(*Step size*)
(*End of the coordinates of the boundaries*)


(*The boundary values of V at the boundaries*)
vx[iini]=-1.0;vx[ilast]=1.0;vy[jini]=-1.0;vy[jlast]=1.0;
(*End of the boundary values of V at the boundaries*)

(*Defining the values of coordinates of x[i],y[j]*)
For[j=jini,j<=jlast,j++,
	y[j]=y[jini]+(j-jini)*Deltax;
];(*end for j*)

For[i=iini,i<=ilast,i++,
	x[i]=x[iini]+(i-iini)*Deltax;
];(*end for j*)
(*End of defining the values of coordinates of x[i],y[j]*)

(*Emptying all values in V[i,j]*)
For[j=jini,j<=jlast,j++,
For[i=iini,i<=jlast,i++,
V[i,j]=0;
Vold[i,j]=0;
];(*end for i*)
];(*end for j*)
(*End of emptying all values in V[i,j]*)


(*Overwritting V at the boundaries with the boundary conditions*)
For[j=jini,j<=jlast,j++,
For[i=iini,i<=ilast,i++,
If[i==iini,V[i,j]=vx[iini]];(*end if i=iini*)
If[i==ilast,V[i,j]=vx[ilast]];(*end if i=ilast*)
If[j==jini,V[i,j]=vy[jini]+((vy[ilast]-vy[iini])/(y[ilast]-y[iini]))*(i-iini)*Deltax];
(*end if j=jini*)
If[j==jlast,V[i,j]=vy[jini]+((vy[ilast]-vy[iini])/(y[ilast]-y[iini]))*(i-iini)*Deltax];
(*end if j=jlast*)
];(*end for i*)
];(*end for j*)
(*Overwritting V at the boundaries with the boundary conditions*)


tV01=Table[{x[i],y[j],V[i,j]},{i,iini,ilast},{j,jini,jlast}];
tV02=Flatten[tV01,1];
Print["Print from initialisation.m"];
Print[ListContourPlot[tV02,Mesh->None,ContourLabels->All]];

(*Check that indeed the initial guess generated by the code is according to what we want*)
For[j=jini,j<=jlast,j++,
Print[Table[V[i,j],{i,iini,ilast}]]
];(*end for j*)
