For existence theorem to be used f(x,y) needs to be continuous in the interval
\begin{equation}
R=\left\{(x, y):\left|x-x_{0}\right| \leq a,\left|y-y_{0}\right| \leq b\right\}, \quad(a, b>0)
\end{equation}
To find the interval actually you should randomly pick the area by yourself because you are analysing the equation and you pick those values of a and b according to the function so, you should pick some interval and test it with the existence and uniqueness theorems but they are not fully telling you the interval of validity so you should find the interval first according to some techniques.
The order is always, first one fixes $a,b$, then computes the maximum $M=M(a,b)$ (and confirms the Lipschitz condition), then derives $h=\min(a,\frac{b}{M(a,b)})$ as the interval radius where the solutions can not leave the box/cylinder via its horizontal boundaries.
Now it may happen that the initial length parameters $a,b$ where chosen so fortunate that $h=a=\frac{b}{M(a,b)}$ by this mechanism. If $M(a,b)$ is an easy enough expression, one can help this luck by some computation. As there are two free parameters, one can then optimize $b$ so that $a$ becomes as large as possible.
But none of that is necessary for the theoretical application of this construction in the existence theorem, and usually the results of the above optimization still fall (very) short of other estimates of intervals contained in the maximal domain.
I do not see the difference between the two methods you describe. The first is perhaps a little more formal. Note that if you do this step-by-step, then $M(h,b)$ will most often be smaller than $M(a,b)$, which means that one could relax $h$ a little in direction $a$ and restart with a new $a$ value, etc. Implicitly this is a fixed-point equation for $a$ where the value of $b$ remains unchanged. Having a solution or sufficient approximation of this giving $a=g(b)$ as a function of $b$ with $g(b)\le\frac{b}{M(g(b),b)}$, one can then proceed to find an optimal $b$ as a maximum of $g$.
In your written example calculations, there is no difference between the methods in the cases of the autonomous equations, as $a$ does not influence $M$, so that $g(b)=\frac{b}{M(\cdot,b)}$ directly. In the other cases, I think that the examples are simple enought that the admissible set $\{(a,b):aM(a,b)\le b\}$ is convex so that it does not matter if you approach the boundary in horizontal direction like in method 2 or in vertical direction like in method 1. I do not have counter-examples, but there might be situations where this convexity fails and you get different results from both methods.
Best Answer
In your case $x_0 = y_0 = 0$. For any $a > 0, b > 0$, the function $f(x, y) = y^2 + \cos^2 x$ is defined and Lipschitz continous on the rectangle $R=\{|x|\leq a, |y|\leq b\}$.
On this rectangle, $$|f(x, y)| \le |y|^2 + |\cos x|^2 \le b^2 + 1$$ with equality for $x = 0$ and $y = b$, so the supremum is $M = b^2 + 1$ and therefore $$h= \min\{a,\frac{b}{b^2 + 1}\}$$
The Picard existence theorem states that the IVP has a (unique) solution on the interval $[-h, h]$, or – if you restrict the problem to $x \ge 0$ – on $[0, h]$.
The task is now to choose $a$ and $b$ such that $h$ becomes as large as possible. From the AM-GM inequality is follows that $$ b = \sqrt{1 \cdot b^2} \le \frac{b^2 + 1}2 \Longrightarrow \frac{b}{b^2 + 1} \le \frac 12 $$ with equality for $b=1$.
It follows that $h \le \frac 12$ for any choice of $a, b$, and $h = \frac 12$ for $a = \frac 12, b = 1$. So $h = \frac 12$ is the largest value that can be obtained by this method.