No, there is no better way in the general case.
You can create a single function that calculates the 10000 different values, and you can ask fsolve() to work on it, passing in a vector of 10000 starting values. fsolve() would see this as-if it were a problem with 10000 inter-related variables, and it would spend time building the 10000 by 10000 gradient at each step, trying to figure out how to solve all of the equations as-if each equation depended on every variable. This might produce solutions, but it is not going to produce solutions efficiently.
In some cases you can use calculus to analyze the equation with regards to variable c and get hints about where the solution must be based upon the value of c. For multinomials (polynomials extended to multiple variables) this can be very effective. For example, the family fun(x,c) = x^2 + c*x + 1 easily reduces down to a pair of roots in terms of c, and then you can substitute in the list of c values to get the list of roots without having to solve each one individually.
