You can directly compute the optimal strategy and its expected value with a dynamic program (backward induction).
Consider the possible states of the game, which can be completely described by the (number of -1 cards remaining, number of +1 cards remaining), with 16 possibilities.
Arrange them in a square grid as follows (it might be better on paper if you make it a diamond with (3,3) on the left and (0,0) on the right)
$$\begin{array}{ccccccc}
\stackrel{(0)}{(3,3)} & \rightarrow & \stackrel{(1)}{(3,2)} & \rightarrow & \stackrel{(2)}{(3,1)} & \rightarrow & \stackrel{(3)}{(3,0)}\\
\downarrow && \downarrow && \downarrow && \downarrow\\
\stackrel{(-1)}{(2,3)} & \rightarrow & \stackrel{(0)}{(2,2)} & \rightarrow & \stackrel{(1)}{(2,1)} & \rightarrow & \stackrel{(2)}{(2,0)}\\
\downarrow && \downarrow && \downarrow && \downarrow\\
\stackrel{(-2)}{(1,3)} & \rightarrow & \stackrel{(-1)}{(1,2)} & \rightarrow & \stackrel{(0)}{(1,1)} & \rightarrow & \stackrel{(1)}{(1,0)}\\
\downarrow && \downarrow && \downarrow && \downarrow\\
\stackrel{(-3)}{(0,3)} & \rightarrow & \stackrel{(-2)}{(0,2)} & \rightarrow & \stackrel{(-1)}{(0,1)} & \rightarrow & \stackrel{(0)}{(0,0)}\\
\end{array}
$$
The entries are on top (score if you stop here), bottom (#-1, #+1) remaining in the deck.
The trick is to work backwards from the (0,0) corner and decide at each state whether you want to continue or not. Examples:
- There is no decision at (0,0), it is worth 0.
- At (0,1) the choice is between taking -1, or drawing a card which gets 0. Since drawing is better, we now know (0,1) is also worth 0.
- At (1,0) we take 1 instead of drawing which gets 0.
- At (1,1) a real decision comes up. Stopping is worth 0. Drawing gets you 1/2 chance to move to (1,0) [worth 1] and 1/2 chance to move to (0,1) [worth 0]. So drawing is worth 1/2 on average, and it is optimal to do so.
You can continue filling in all the states to find the optimal strategy. Note that unequal card counts matter: at say (1,2), drawing gives you 1/3 chance to move to (0,2) and 2/3 chance to move to (1,1).
The filled in square looks like:
$$\begin{array}{ccccccc}
\stackrel{17/20}{(3,3)} & \rightarrow & \stackrel{6/5}{(3,2)} & \rightarrow & \stackrel{\mathbf{2}}{(3,1)} & & \stackrel{\mathbf{3}}{(3,0)}\\
\downarrow && \downarrow && &&\\
\stackrel{1/2}{(2,3)} & \rightarrow & \stackrel{2/3}{(2,2)} & \rightarrow & \stackrel{\mathbf{1}}{(2,1)} & \rightarrow^{?} & \stackrel{\mathbf{2}}{(2,0)}\\
\downarrow && \downarrow && \downarrow^{?} &&\\
\stackrel{1/4}{(1,3)} & \rightarrow & \stackrel{1/3}{(1,2)} & \rightarrow & \stackrel{1/2}{(1,1)} & \rightarrow & \stackrel{\mathbf{1}}{(1,0)}\\
\downarrow && \downarrow && \downarrow &&\\
\stackrel{0}{(0,3)} & \rightarrow & \stackrel{0}{(0,2)} & \rightarrow & \stackrel{0}{(0,1)} & \rightarrow & \stackrel{\mathbf{0}}{(0,0)}\\
\end{array}$$
States where you stop have their value bolded. At (2,1) it doesn't matter if you draw or stop.
Since you have made value-maximizing choices at every step including the effects of later choices, Strategy 2 is proven optimal, with value exactly 17/20.
Letting $E_N$ be the optimal expected value for even $N$, we can show that
$$
E_N\sim c\sqrt N
$$
for $N$ large, where $c\approx 0.369136$.
I will give a rough argument, which it should be possible to make rigorous by going through the details. Define a process $X^N_t$ for $0\le t\le 1$ as follows. Let $X^N_{k/N}$ be the number of black cards minus the number of red cards after $k$ have been selected. We interpolate in the range $[(k-1)/N,k/N)$ as constant. This is a random walk conditioned on $X_1=0$. Scaling, $N^{-1/2}X^N_t$ converges, in distribution, to a limit $X$ which is a standard Brownian bridge on the interval $[0,1]$ (i.e., a standard Brownian motion conditioned on $X_1=0$). So, commuting the limit with the supremum,
$$
N^{-1/2}E_N=\sup_\tau E[N^{-1/2}X^N_\tau]\to\sup_\tau E[X_\tau].
$$
The supremum is taken over stopping times $0\le\tau\le 1$.
So, $E_N\sim c\sqrt N$ where $c=\sup_\tau E[X_\tau]$ is the solution to the continuous-time version of the game described.
To solve for $c$, let $f(t,x)$ be the solution starting from $X=x$ at time $t$. That is,
$$
f(t,x)=\sup_{\tau\ge t} E[X_\tau\vert X_t=x].
$$
By considering $\tau=t$ and $\tau=1$, we see that $f(t,x)\ge\max(x,0)$. We need to compute $c=f(0,0)$. We can reduce the dimensionality of the problem: by a simple symmetry, for each fixed time $t$, the process $Y_s=(1-t)^{-1/2}X_{t+s(1-t)}$ is also a Brownian motion conditioned on $Y_1=0$. So,
$$
\begin{align}
f(t,x)&=\sup_{\tau} E[X_{t+\tau(1-t)}\vert X_t=x]\\
&=\sup_\tau(1-t)^{1/2} E[Y_\tau\vert Y_0=(1-t)^{-1/2}x]\\
&=(1-t)^{1/2}g((1-t)^{-1/2}x)
\end{align}
$$
where I am setting $g(x)=f(0,x)$. Next, the Brownian bridge satisfies the Stochastic differential equation
$$
dX_t = dW_t - (1-t)^{-1}X_tdt
$$
for standard Brownian motion $W$.
On the domain where $f(t,x) > x$ then $f(t,X_t)$ is locally a martingale, so it satisfies the Kolmogorov backward equation (PDE),
$$
\frac{\partial}{\partial t}f+\frac12\frac{\partial^2}{\partial x^2}f-(1-t)^{-1}x\frac{\partial}{\partial x}f=0.
$$
Substituting in the expression for $f$ in terms of $g$ gives a linear ODE
$$
g^{\prime\prime}(y)-yg^{\prime}(y)-g(y)=0.
$$
This has the general solution
$$
g(y)=e^{y^2/2}(a+bN(y))
$$
for constants $a,b$, where $N(y)=(2\pi)^{-1/2}\int_{-\infty}^ye^{-t^2/2}dt$ is the cumulative normal distribution function. As $g$ is increasing, it cannot blow up in the limit $y\to-\infty$, so we have $g(y)=be^{y^2/2}N(y)$ on the domain $g(y) > y$. Plugging in $y=0$ gives $c=g(0)=b/2$, so $g(y)=2ce^{y^2/2}N(y)$ over $g(y) > y$. As $g(y)$ cannot blow up faster than linearly in $y$, there must be a $y_0 > 0$ for which $2ce^{y_0^2/2}N(y_0)=y_0$, so that $g(y)=y$ for $y > y_0$. The optimal value of $c$ for which this holds is given by
$$
c=\sup_{y > 0}\frac{ye^{-y^2/2}}{2N(y)}.
$$
By differentiating, the function of $y$ on the right has a unique local maximum. I am not sure if there is an analytic solution, but I tried in wolfram alpha. Expressing in terms of the error function,
$$
c=\sup_{y > 0}\frac{ye^{-y^2/2}}{1+{\rm erf}(y/\sqrt2)}
$$
we get $c\approx0.369136$.
For comparison, I also computed the $\alpha$ in the linked paper as $\alpha\approx0.839923675692373$, so $c=(1-\alpha^2)(\pi/2)^{1/2}\approx 0.369136$ and we are in agreement.
Best Answer
You should never stop when you are losing because you can guarantee $0$ by drawing all the cards. Clearly you should stop after three cards if you are $+1$ or after two cards if you are $+2$ as you can only get worse. You should not stop if you are even after two because you can only get better. The only question is whether to draw if you are $+1$ on the first draw.
We compute the expectation if you draw red first and draw again. You have $\frac 13$ chance of drawing red again and ending $+2, \frac 13$ chance of drawing two blacks next and ending $0$, and $\frac 13$ chance of drawing black-red and ending $+1$. This gives $+1$, so it doesn't matter whether you draw or not.
Now we compute the expectation at the start. If you draw red (probability $\frac 12$) you end $+1$. If you draw black and then draw two reds (probability $\frac 16$) you end $+1$ Otherwise you break even with probability $\frac 13$. Overall, the value is $\frac 23$