You should realize that there is an intrinsic limitation of sparse-direct methods like this for matrices that come from large 3d meshes/grids — they scale poorly for 3d problems, and tend to run out of memory before the problems become large. You could try a different sparse-direct solver like MUMPS (GitHub - JuliaSmoothOptimizers/MUMPS.jl: A Julia Interface to MUMPS or GitHub - JuliaInv/MUMPSjInv.jl: MUMPS interface for Julia), which might suffice in the short term, but in the long term you will hit the same scaling issues.
Here is a nice slide on the scaling (by Per Persson) that I use in one of my classes:
Note that for an n \times n \times n grid, N=n^3 so the storage scales as O(n^4). For more information, see e.g. Tim Davis’ book.
In consequence, for large 3d PDEs, the sparse-matrix problems must ultimately be solved by iterative methods (e.g. GMRES, CG, etc.) rather than sparse-direct methods. There are a few packages for iterative solvers in Julia, e.g. IterativeSolvers.jl, but to use iterative solvers effectively you have to know a bit more, and you often want to find a problem-specific preconditioner.