num=64(* Rozmiar siatki *)

n=(num-2)^2 (* Rozmiar macierzy *)

generatorMacierzy[num_]:=Module[{diagN},
diagN=ConstantArray[Join[Table[-1,num-3],{0}],num-2]//Flatten//Most;
SparseArray[{
Band[{1,1}]->4,
Band[{1,2}]->diagN,
Band[{2,1}]->diagN,
Band[{1,num-1}]->-1,
Band[{num-1,1}]->-1},{(num-2)^2,(num-2)^2}]
];

generatorRHS[num_,f_]:=
Module[{F,i,j},
F=Table[ToExpression["fx"<>ToString[i]<>"y"<>ToString[j]],{i,0,num-1},{j,0,num-1}];
Table[F[[i,1]]=f[Pi(i-1)/(num-1),0],{i,2,num-1}];
Table[F[[i,num]]=f[Pi(i-1)/(num-1),Pi],{i,2,num-1}];
Table[F[[1,j]]=f[0,Pi(j-1)/(num-1)],{j,1,num}];
Table[F[[num,j]]=f[Pi,Pi(j-1)/(num-1)],{j,1,num}];
Table[
If[i==2&&j==2,F[[1,2]]+F[[2,1]],
If[i==2&&j==num-1,F[[1,num-1]]+F[[2,num]],
If[i==num-1&&j==num-1,F[[num-1,num]]+F[[num,num-1]],
If[i==num&&j==2,F[[num,2]]+F[[num-1,1]],
If[i==2,F[[1,j]],
If[j==2,F[[i,1]],
If[j==num-1,F[[i,num]],
If[i==num-1,F[[num,j]],
0
]
]
]
]
]
]
]
]
,{i,2,num-1},{j,2,num-1}]//Flatten
];

m=generatorMacierzy[num] (* Macierz 9x9 *)

(*MatrixForm[m]*)

Normal[m] (* To można przekleić do kodu w C *)

f=Function[{x,y},Cosh[y]Sin[x]];

generatorRHS[num,f];

b=N[%,16]

Import["rozwLAPACK64.m"]

LinearSolve[m,b]

rozwLAPACK

%-%%

%/$MachineEpsilon

Abs[%]//Max

Table[f[Pi/(num-1)*(i-1),Pi/(num-1)*(j-1)],{i,2,num-1},{j,2,num-1}]//Flatten//N;

LinearSolve[m,b];

%-%%

Total[Abs[%]]

Mean[Abs[%%]]