I am working on an FEM program written in C, for my undergraduate degree which needs very large arrays (arrays consisting a[1 000 000] elements) to store data, and then manipulating them. It uses 2D arrays also, which have similar unusually huge sizes (say a[100 000][100 000]).
The program terminates without showing any error. It has been found that just when the program will execute the routine for generating mesh (which uses large 2D arrays), it crashes.
When you reach such large sizes you should be thinking whether or not your matrix really has 100 000 x 100 000 elements or if most of them are null. If most are zeros you should use sparse matrices. That should alleviate memory use.
After that you should try and use matrix decomposition (such as lower upper) to solve your system, I believe you should be able to find implementations in your favourite language.
There are other people interested in such large systems, so try and see how they did it, and take advantage of approximation / iterative solvers.