Maximizing a function by finding derivative I want to find the value of $\vec{p}$, $p_s$, $p_t$ each of which is a function of the form $f:\mathbb{R}^2 \to \mathbb{R}$ that maximize the following function :
$$\begin{align} 
\int_\mathbb{R^2} \{ p_s &- \alpha(x,y)|\vec{p}| - \beta(x,y)(p_s - C_s) \\
&- \gamma(x,y)(p_t - C_t) + \lambda(x,y)(\vec{\nabla}\cdot\vec{p} - p_s + p_t)\}
\end{align} $$
Each of $\alpha, \beta, \gamma$ is constrained to be non negative, $\lambda$ is unconstrained. $C_s$ and $C_t$ are given functions of the form $f:\mathbb{R}^2 \to \mathbb{R}$. 
Edit : 
Also we want the sup to be finite and therefore can put constraints on $\alpha, \beta, \gamma, \lambda$ accordingly. 
The problem is a convex optimization problem. 
I am hoping to minimize it by using a derivative approach : minimize function for a given variable by finding a derivative with respect to the variable and equating it to zero. I just dont know how to do that here. 
 A: There is no maximum unless $\beta+\lambda\equiv1$ and $\gamma\equiv\lambda$. The terms containing $p_s$ add up to
$$\iint_{\mathbb R_2}p_s(1-\beta-\lambda)\mathrm dx\mathrm dy\;,$$
so you can make the integral arbitrarily large by choosing $p_s$ arbitrarily large (positive or negative depending on the sign of $1-\beta-\lambda$). Likewise, the terms containing $p_t$ add up to
$$\iint_{\mathbb R_2}p_t(\lambda-\gamma)\mathrm dx\mathrm dy\;.$$
This fixes $\beta$ and $\gamma$ in terms of $\lambda$, and then $p_s$ and $p_t$ no longer occur in the integral. We can drop the constant terms containing $\beta$ and $\gamma$, leaving us with
$$I=\iint_{\mathbb R^2}\left[- \alpha|\vec{p}|+ \lambda(\vec{\nabla}\cdot\vec{p})\right]\mathrm dx\mathrm dy\;.$$
Let's restrict ourselves to solutions $\vec{p}$ that decay sufficiently rapidly at infinity that we can integrate by parts and omit the boundary term:
$$I=\iint_{\mathbb R^2}\left[- \alpha|\vec{p}|- (\vec{\nabla}\lambda)\cdot\vec{p}\right]\mathrm dx\mathrm dy\;.$$
So we want to find $\vec{p}$ such that its magnitude (weighted by $\alpha$) is small but its component along $-\vec{\nabla}\lambda$ is large. This will work best if we choose $\vec{p}$ to always point along $-\vec{\nabla}\lambda$, i.e. $\vec{p}=-\mu\vec{\nabla}\lambda$, with as yet undetermined $\mu(x,y)$. With $\vec{g}:=\vec{\nabla}\lambda$, this gives
$$
\begin{eqnarray}
I
&=&
\iint_{\mathbb R^2}\left[- \alpha|-\mu\vec{g}|- \vec{g}\cdot(-\mu\vec{g})\right]\mathrm dx\mathrm dy
\\
&=&
\iint_{\mathbb R^2}\left[- \alpha|\mu||\vec{g}|+\mu|\vec{g}|^2\right]\mathrm dx\mathrm dy
\\
&=&
\iint_{\mathbb R^2}|\vec{g}|\left[- \alpha|\mu|+\mu|\vec{g}|\right]\mathrm dx\mathrm dy
\;.
\end{eqnarray}$$
Thus, we must have $\alpha\ge|\vec{\nabla}\lambda|$ everywhere; else we could make the integral arbitrarily large by choosing $\vec{p}$ to vanish wherever that inequality holds and to make the integral arbitrarily large where it doesn't.
Given that inequality, the integral is forced to be non-positive, and thus it is maximized for $\vec{p}\equiv0$.
The first conclusion, $\alpha\ge|\vec{\nabla}\lambda|$, holds independent of the assumption about the behaviour at infinity, since we can choose $\vec{p}$ to have that behaviour if we want. The second conclusion, $\vec{p}\equiv0$, however, holds only under that assumption, since the boundary term might compensate a negative value of the integral. I'm not sure how to treat the general case, without decay at infinity.
[Edit:] Quite generally speaking, the parts of the integral that depend on $\vec{p}$ all scale linearly if you scale $\vec{p}$ by a positive scale factor, so there can't be any maximum other than $\vec{p}=0$ -- all you can do by imposing constraints on $\alpha$, $\beta$, $\gamma$ and $\lambda$ is to ensure that $\vec{p}=0$ is indeed a maximum.
