Why are recursives methods useful for FIR filter design?
Sure, FIR filters are "simple" linear objects. But a typical FIR filter design task looks like this:


  Find a set of numbers (the filter coefficients) such that the
  magnitude of their Discrete Fourier Transform is as close as possible to a
  target response. "As close as possible" being defined as "by
  minimizing the $L_\infty$ norm of the difference between ideal and actual
  response over some frequency ranges".


This is a non-linear, sometimes non-convex optimization problem and there is no universal, trivial procedure to solve those.

We have to use procedures like the Parks-McClellan algorithm because the real-world engineering constraints we want to impose when designing FIR filters do not translate into simple linear or quadratic constraints. For example, a least-square error criterion (L2 norm often leads to more tractable procedures than other norms...) would not be an acceptable choice for engineering applications, because it could allow solutions that are close to the required response on average, but with an outlier value.