# this is the boundary control example with alpha = 0
#
param n default 39;
param n1 := n+1;
param a default .0;
param g default 0;
param t default 5;

var x{0..n1, 0..n1};

param h := 1/n1;
param h2 := h^2;

param z{i in 1..n, j in 1..n} := 1;

minimize f:	.5*h2*sum{i in n1/4..3*n1/4, j in n1/4..3*n1/4} 
(x[i,j]-z[i,j])^2 + .5*h*a*sum{i in 1..n} x[i,n1]^2;

s.t. pde{i in 1..n, j in 1..n}:
	4*x[i,j] - sum{k in {-1,1}} (x[i+k,j] + x[i,j+k]) = g*h2;

s.t. bc1 {i in 1..n}: x[i,0] = x[i,1];
s.t. bc2 {i in 1..n}: x[0,i] - x[1,i] = h * (x[0,i] - t);
s.t. bc3 {i in 1..n}: x[n1,i] - x[n,i] = h * (x[n1,i] - t) ;

s.t. sc{i in 0..n1, j in 0..n}: 0<= x[i,j] <= 
                      if n1/4<=i<=3*n1/4 && n1/4<=j<=3*n1/4 then 3.15 else 10;

s.t. cc{i in 1..n}: 0 <= x[i,n1] <= 10;    # control constraints

s.t. fix_corners1{j in {0}, i in {0,n1}}: x[i,j] = 0;
s.t. fix_corners2{j in {n1}, i in {0,n1}}: x[i,j] = 10;