I am trying to maximize this function on the "p" value:

```
g2ces(p_, theta_, I_, a_, c1_) :=
theta*p*Min(I,
a*p^c1) + (theta*
theta)*(((I - Min(I, a*p^c1))/a)^(1/c1)*(I -
Min(I, a*p^c1)) + ((I/a)^(1/c1))*I)
```

In the following values of other parameters (which I take as an example), we obtain:

```
FindMaximum({g2ces(p, 0.5, 2, 2.5, -2)}, p)
```

p = 1.25. This is the right solution.

And plot this on several "p" values gives:

```
Plot(Evaluate@
Table(g2ces(p, theta, 2, 2.5, -2), {theta, 0.5, 0.5,1}), {p, 0, 5})
```

This solution is derived analytically and is a special case. This function is equivalent to the following function, "VDPces", which is easier to generalize (and, therefore, I want to work with this version of the function instead of the "g2ces"):

```
Rev(p_, I_, a_, c1_) := p*Min(I, a*p^c1)
V2ces(I_, a_, c1_) := MaxValue({Rev(p, I, a, c1)}, p)
VDPces(p_, theta_, I_, a_, c1_) :=
theta*(Rev(p, I, a, c1) +
theta*V2ces(I - Min(I, a*p^c1), a, c1)) + (1 - theta)*(0 +
theta*V2ces(I, a, c1))
```

When I maximize this function, I get a different solution:

```
FindMaximum({VDPces(p, 0.5, 2, 2.5, -2)}, p)
```

p = 1,111803 (incorrect)

and a different plot

```
Plot(Evaluate@
Table(VDPces(p, theta, 2, 2.5, -2), {theta, 0.5, 0.5, 1}), {p, 0,
5})
```

That seems like:

Obviously something goes wrong with the optimization.