Question: I implemented a finite differences algorithm to solve a PDE.

The grid is a structured 2D domain of size [Nx, Nz], solved Nt times.

I pre-allocate the object containing all solutions:

```
sol = zeros(Nx, Nz, Nt, 'single') ;

```

This becomes very easily too large and I get a '*out of memory*' error.
Unfortunately `sparse` doesn't work for N-dimensional arrays.

For the sake of the question it's not important to know the values, it goes without saying that the RAM usage grows exponentially with decreasing the grid spacing and increasing the simulation time.

I am aware that I do **not** need to store each time instant for the purpose of the advancement of the solution. It would be sufficient to just store the previous two time steps. **However**, for post-processing reasons I need to access the solution at all time-steps (or at least at a submultiple of the total number).It might help to specify that, even after the solution, the grid remains predominantly populated by zeros.

Am I fighting a lost battle or is there a more efficient way to proceed (other type of objects, vectorization...)?

Thank you.

Answer: 
You can create a a cell array of sparse matrices to store the results. However computations can be performed on full matrices if working with a full matrix is faster than sparse matrix and convert the full matrix to sparse matrix and place it in the cell.
You could store the array in sparse, **linear** form; that is, a column vector with length equal to the product of dimensions:

```
sol = sparse([], [], [], Nx*Nz*Nt, 1); % sparse column vector containing zeros

```

Then, instead of indexing normally,

```
sol(x, z, t),

```

you need to translate the indices `x`, `z`, `t` into the corresponding linear index:

* For **scalar** indices you use

```
sol(x + Nx*(z-1) + Nx*Nz*(t-1))

```

You can define a helper function for convenience:

```
ind = @(sol, x, y, t) sol(x + Nx*(z-1) + Nx*Nz*(t-1))

```

so the indexing becomes more readable:

```
ind(sol, x, z, t)

```
* For **general** (array) indices you need to [`reshape`](https://www.mathworks.com/help/matlab/ref/reshape.html) the indices along different dimensions so that [implicit expansion](https://blogs.mathworks.com/loren/2016/10/24/matlab-arithmetic-expands-in-r2016b/) produces the appropriate linear index:

```
sol(reshape(x,[],1,1) + Nx*(reshape(z,1,[],1)-1) + Nx*Nz*(reshape(t,1,1,[])-1))

```

which of course could also be encapsulated into a function.

**Check** that the conversion to linear indexing works (general case, using non-sparse array to compare with normal indexing):

```
Nx = 15; Nz = 18; Nt = 11;
sol = randi(9, Nx, Nz, Nt);
x = [5 6; 7 8]; z = 7; t = [4 9 1];
isequal(sol(x, z, t), ...
    sol(reshape(x,[],1,1) + Nx*(reshape(z,1,[],1)-1) + Nx*Nz*(reshape(t,1,1,[])-1)))

```

gives

```
ans =
  logical
   1

```