I have completed part(a) to find the $y$-Lipschitz constant $L=3$, I'm just not sure how to start part (b), determining a step size to guarantee a global error less than $10^{-3}$.
Any help would be greatly appreciated.
Thanks!
(after comments) I have followed the link and put my values into the formula $$|w_i-y_i|\le \frac{Mh}{2L}(e^{L(x_i-a)}-1)$$ with $w_i−y_i=0.001$, $M=40$, $L=3$ and $x_i−a=3−1=2$; but then when everything cancels out, I get $h=3.727\cdot 10^{-7}$. Any idea where I might be going wrong?
Solving the problem with 5 and 10 steps gives the table
where $y_1=y_{2h}(x)$ are the values for $5$ steps with $2h=0.4$, $y_2=y_h(x)$ the corresponding values for a solution with $h=0.2$ and $c(x)=(y_{2h}(x)-y_h(x))/h$ is the estimate for the leading coefficient in $$y_h(x)=y_{\rm exact}(x)+c(x)h+O(h^2).$$ So with $|c(x)|\le 3$ one would need $h=10^{-3}/3$ for the goal of the task.
The theoretical bound will lead to a radically smaller bound for $h$.
(after the discussion in comments below the question)
The value you get from the cited formula is correct, I got the same value for my last remark above. Note that the theoretical bound is in all steps of its derivation strictly pessimistic to apply to all possible situations. In practice the contributions are not always maximum and some local errors may have opposite sign and thus partially cancel out.
With some care one can show that $|y′′(x)|≤12$, but that improves the bound only by a factor of about $3$. At all points $f(x,y)\in[1,3]$, so for any solution $x-1\le y(x)-y_0\le 3(x-1)$, next $y''(x)=Df(x,y)=f_x+f_yf=-\cos(xy)(y+xf(x,y))$, so on this wedge you get $Df\in [-1,1]\cdot([-0.5,5.5]+[1,2]\cdot[1,3])\subset[-11.5,11.5]$.
My answer contains an a-posteriory estimate of the error behavior, essentially extrapolated from 2 points, but still surprisingly accurate. Plotting the estimates for a range of step sizes shows how stable this estimate is and how it evolves with $x$. As one can see, there is no visible exponential growth in this span of the solution.