That might work, but it has far more variables than equations. I recommend thinking instead as follows:
The region $ 0 \leq x, y \leq 1$ is a square. More importantly, it's a compact region. So we know that our function will take a maximum and a minimum somewhere. So it will either attain its extrema in the interior, or on the boundary.
To check the interior, perform the usual max/min test: take the gradient and set to zero. Either use the multivariable 'second-derivative test' or just compare the points to check for maxima and minima. Now, you should check the boundaries. This means you could do the regular Lagrange multipliers method 4 times, one with each constraint
$$\begin {align} y &= 0; \quad x = 0 \\ y &= 0; \quad x = 1 \\ y &= 1; \quad x = 0 \\ y &= 1; \quad x = 1 \end{align}$$
I want to emphasize that I would do these constraints separately rather than together. Each one is very trivial to solve - but you only pay attention to solutions within your region of interest. There is one more place that you need to check, and that's the corners (here, the boundary of the boundary: it's as if we treated our domain as 4 separate lines, themselves closed and bounded).
Finally, you compare these three areas: the interior, the line-boundaries, and the corner boundaries. Because the region is compact, there will be an absolute max and absolute min (perhaps multiple). That seems to be the easiest way to do this, especially as all the equations in this case are very simple to solve (all linear).
I should note that you do not need to use Lagrange multipliers if you don't want. You could interpret the lines as paths, and maximize/minimize the curve along each of the paths (which does not require lagrange, but instead a parameterization). Then you still compare the interior and the boundary. But this sounds not so fun to me, and not worth it in this case.
What goes wrong is that the minimum might happen to be at a point where a derivative might not exist. Try e.g. minimize $y$ subject to $g(x,y) = y - |x| = 0$. The minimum is at $(0,0)$, where $\frac{\partial g}{\partial x}$ does not exist.
Best Answer
Here is a solution which is a mixture of elimination and Lagrange. Using the only constraint that's an equality, we can substitute $z = 1-3x-2y$ into the function and the other constraints. So now the problem is to maximize $$ f(x,y) = (x+y+(1-3x-2y))^3 = (1-2x-y)^3 $$ subject to the constraints $x\geq 0$ and $$ x^2 + y^2 \leq 1-3x-2y $$ The first inequality describes the right half of the $xy$-plane. Complete the square in the second inequality: \begin{align*} x^2 + 3x + y^2 + 2y &\leq 1 \\ x^2 + 3x + \frac{9}{4} + y^2 + 2y + 1 &\leq \frac{17}{4} \\ \left(x+\frac{3}{2}\right)^2 + \left(y + 1\right)^2 &\leq \frac{17}{4} \end{align*}
Therefore, this inequality describes a disk of radius $\frac{\sqrt{17}}{2}$ centered at $\left(-\frac{3}{2},-1\right)$. The two together form the constraint set. It's the part that's shaded red and blue below.
When looking for the extreme points, you want to search for ordinary critical points inside the constraint set, and constrained critical points on the boundary of the set. For the first, we use the first-order conditions.
\begin{align*} \frac{\partial f}{\partial x} &= -6(1-2x-y)^2 \\ \frac{\partial f}{\partial y} &= -3(1-2x-y)^2 \\ \end{align*} The only way both of these can be zero is if $1-2x-y=0$, or $2x+y = 1$. This describes a line through $(0,1)$ and $(1/2,0)$, shown in green on the diagram below. You can see it doesn't touch the constraint set at all. So there are no unconstrained critical points to check.
The boundary of the constraint set includes two curves $x=0$ and $\left(x+\frac{3}{2}\right)^2 + \left(y + 1\right)^2 = \frac{17}{4}$. First, use Lagrange to find any critical points along $x=0$. We get \begin{align*} -6(1-2x-y)^2 &= \lambda \\ -3(1-2x-y)^2 &= 0 \\ x&= 0 \end{align*} The third and second equation tell us that $(1-y)^2 =0$, or $y=1$. Again, this is outside the constraint set. Next, use Lagrange with the second constraint equation: \begin{align*} -6(1-2x-y)^2 &= 2\lambda \left(x+\frac{3}{2}\right) \\ -3(1-2x-y)^2 &= 2\lambda \left(y + 1\right) \\ \left(x+\frac{3}{2}\right)^2 + \left(y + 1\right)^2 &= \frac{17}{4} \end{align*} The first two equations tell us that $\lambda (2x+3) = \lambda(4y+4)$, so either $\lambda = 0$ (this forces $1-2x-y=0$, which has already been ruled out) or \begin{align*} 2x+3 &= 4y+4 \\ y &= \frac{1}{2}x-\frac{1}{4} \end{align*} This line, shown in orange below, does intersect the boundary of the constraint set in two places: $\left(0,-\frac{1}{4}\right)$ (when $x=0$) and $\left(\sqrt{\frac{17}{5}}-\frac{3}{2},\frac{1}{2}\sqrt{\frac{17}{5}} - 1\right)$.
Finally, because this constraint set is not smooth but has corners, we have to include those in our check points. They occur when $x=0$ and $\left(x+\frac{3}{2}\right)^2 + \left(y + 1\right)^2 = \frac{17}{4}$, or $\left(0,-1 \pm \sqrt{2}\right)$.
So, four points to check:
The greatest of these values by far is the one at $\left(0,-1-\sqrt{2}\right)$. This is pretty messy, so I'm not 100% confident in the calculations, but Wolfram Alpha agrees on the maximum.