Operator learning has emerged as a transformative paradigm for modeling complex dynamical systems, yet traditional neural operators face challenges in computational efficiency, scalability, and integration with system-level numerical analysis. This talk introduces RandONets - a novel class of randomized neural operators - that builds on a newly established universal approximation theorem for nonlinear operators with randomized, fixed-weight features. RandONets employ shallow architectures with randomized projections (e.g., Johnson-Lindenstrauss, Random Fourier Features) and a decoupled branch-trunk training to enable fast, accurate, and interpretable operator learning. By reformulating training as a convex least-squares problem, RandONets achieve orders-of-magnitude improvements in training speed and accuracy over DeepONets, as demonstrated on parametric PDE benchmarks approximating evolution and solution operators. Furthermore, Neural Operators have so far primarily been employed as surrogate models to explore the dynamical behaviour through brute-force temporal simulations. Their potential for systematic rigorous system-level analysis-such as bifurcation analysis, remains largely unexplored. We show how to enable RandONets to perform equation-free, system-level analysis of multiscale dynamical systems. Leveraging matrix-free iterative methods in the Krylov subspace, thus enabling neural operators to directly perform tasks and analysis such as steady-state computation, stability analysis and bifurcation tracking - all without recursive temporal simulations. This framework avoids error accumulation and enables direct analysis of high-dimensional systems. We illustrate our framework via three non-linear PDEs: (i) a modified Bratu-Gelfand PDE exhibiting a saddle-node bifurcation; (ii) a one-dimensional Allen-Cahn equation which undergoes multiple concatenated pitchfork bifurcations; and (iii) the FitzHugh-Nagumo (FHN) equation, described by two coupled PDEs that exhibit both a Hopf and a saddle-node bifurcation.