I just took a few minutes to modify the example from this page
http://www.sagemath.org/doc/reference...
I apologize in advance if there is a mistake. Basically there isn't a function to maximize with constraint. So what we can do is minimize the negative of the function. The constraints are of the for '>0' for example c_1 tells the solver that x>0, c_2 tells that -x+1 > 0 => x<1 and so on

```
f = lambda p: - p[0]*p[1] + p[1]*p[1]
c_1 = lambda p: p[0]
c_2 = lambda p: - p[0] + 1
c_3 = lambda p: p[1]
c_4 = lambda p: -p[1] + 2
a = minimize_constrained(f,[c_1, c_2,c_3,c_4],[0.5,0.5])
view(a)
```

The answer I get is (1.0,0.5001)

If I plot the negative of the function in the range you are interested, the minimum looks like it is around that.

Let me know if this works.