I have a procedurally defined Hermitian matrix $M$, i.e. I can get any matrix element by calling a black box function (e.g. a library function), and a vector $Y$. And I have to solve a system of linear equations:
$M\cdot X=Y$.
But $Y$ is such that having $n$ elements, it takes about half of available RAM, another half would be for $X$, so if I try to store $M$, it'll take $n^2$ space, i.e. even if I double the RAM space, this will still place $n2$ matrix columns/rows into swap.
At the same time, all the algorithms which I've found need saving large amounts ($\geq n^2$) of data while solving the system.
So, the question: are there any algorithms to solve systems of linear equations which don't require me to store the matrix, and still are fast enough (maybe not $O(n^3)$, but at least not much slower)?



This looks like a situation where the Kaczmarz method could work. What you do to maintain an approximate solution and then project cyclically onto the hyperplanes which are given by the $k$th equation. More precisely: If you have the $m$the iterate $X^m$ and use the $k$th equation, then you have the next iterate $$X^{m+1} = X^M + \frac{Y_k  a_k^T\cdot X^m}{\a_k\^2}a_i$$ where $a_k$ is the $k$th row of $A$ and $Y_k$ is the $k$the entry of $Y$. Hence, you only need one row of $A$, one entry of $Y$ and the current iterate $X^m$ to perform one iteration, i.e. $2n+1$ space. Also the iteration complexity is very low (it's $\mathcal{O}(n)$) but you usually need a lot of iterations. This methods is widely used in discrete tomography and is also an instance of the "Projection onto convex sets" method. Recently it has been shown by Strohmer and Vershynin that a randomized version of this method has favorable convergence properties (when you pick each column with a probability proportional to its norm). Also "block iterative" versions work, i.e. you take a hyperplane of higher codimension to project on. So, if you have some memory left, you could also take some more rows of $A$ at once... See, e.g. here. 


Being able to get elements of the matrix isn't very useful (particularly if you don't know where the nonzero elements of the matrix are without checking.) Iterative methods can be useful if you have the ability compute matrix vector products $Mx$. You haven't said whether this is possible. It seems quite likely that there's some special structure to your particular problem that would make it possible to simplify this computation. You haven't told us anything about where the system of equations comes from perhaps if you explained this in some detail we could suggest ways to proceed. 

