Numerical method for invers of function

I have a function for a skewed normal distribution, f(x). For this function I like to find the x values such as max(f(x1)) and f(x2) = constant. Since I can’t calculate the derivative analytically I use numerical methods, where the maximum is straight forward with libraries I found. To find x values corresponding to an arbitrary constant I use a more brutish method where I construct arrays of (x, y) values and iterate through to find the closest y value to the constant which gives me the x value.

My question, is there an more elegant way of doing this?

Cheers!

Yes. Use a numerical optimization library like Scipy or Jaxopt. Or if you want to try an analytical solution, try Oryx.