Simplest way to have numerical FOC of a function?

I have the following function for which I need to find the First Order Conditions with respect to one parameter, let’s say \mathbb{u}:

J(u,v,Y,λ) = sum(skipmissing((Y-(u * v')) .^ 2))/2 + (λ/2) * (norm(u)^2 + norm(v)^2)

(This is in the context of an exercise on Collaborative filtering; a machine learning tecnique)

The exercise is very simple, I can do it by hand. But I am wondering which is the best approach, in real life what one would use. Analytical analysis (SymPy) ?
Or automatic differentiation? But in such case I would get a numerical gradient, while I need the result of the equation \partial L / \partial u = 0.

Then you can solve for the result being 0 using a rootfinder algorithm? But if you actually have a closed form I would just use that.