I have the function $2x^2+y$ and one constraint $x-y^2=1$ and want to find maximum value by lagrange multiplier. Intuitively, I see the point $(2,1)$ satisfies $c$ and have value of $f(2,1)=9$.
Using Lagrange function $f(z)-\lambda c(z)$, I get $\lambda^3-4\lambda^2-1$ thus $\lambda=4.06$. However, when plug this into the derivatives of $x$ and $y$ I got different values in particular $x=\frac{4.06}{4}=1.015$ which is much lower than $2$. Am I missing something here?
Note that your function is unbounded. Your constraint gives $x = 1+y^2$, and plugging this into the objective, you get $2\,(1+y^2)^2 + y$. Now, it is easy to see, that this can attain arbitrary large values. This means that your problem has no global maximum, but it may have a local one.
Further, since the gradient of your constraint is $(1, -2\,y)$, the LICQ is satisfied at every feasible point. This means that every local extremum (may it be a maximum or a minimum) can be found by the Lagrange multiplier method.
As you already found, there is only one solution to Lagrange's method. Moreover, it is easy to see that this is a local minimum (consider, e.g., the second derivative of the reduced objective $2\,(1+y^2)^2+y$ which is positive for $y = -0.123$.
Putting everything together, we find that your problem has no local maximum (and no global one).