Suppose I want to solve a system of equations A*x = b, and I have a function A_func such that A_func(x) = A*x. How could I go about solving the system of equations without constructing the matrix A?

I could construct the matrix A if I really wanted to, but I have written enough in terms of A_func that I would like to avoid it if possible.

I saw some similar topics, but they seemed to be focused on the nuances of specific package implementations.

One option is to use an iterative method such as those from Krylov.jl. Have a look at the docs to choose a method based on the properties of A (if any).

SciMLOperators.jl isn’t perfect yet, but it’s pretty close. It has the most integration with the Julia package system by far too, and it’s only growing. So I’d recommend it.

LinearSolve.jl takes AbstractSciMLOperators and will use that with the various Krylov methods. This also means any SciML solver (ODEs, nonlinear solvers, etc.) will also take AbstractSciMLOperators for Jacobian/Hessian definitions.