setup;

% Initial guess
u = ua;

% Solve state equation K*y = beta .* u
y = K \ ( beta .* u );

for i=1:100

	% Solve adjoint equation K*p = y - yd
	p = K \ (y - yd);

	% Compute search direction
	v = ones((n-1)^2,1);

	I = ( lambda*u + beta.*p > 0 );
	v(I) = ua(I);
	I = ( lambda*u + beta.*p < 0 );
	v(I) = ub(I);
	I = ( lambda*u + beta.*p == 0 );
	v(I) = u(I);


	% Compute step lenght s
	w = K \ ( beta .* v);
	s = -( (y - yd)'*M*(w - y) + lambda*u'*M*(v-u) )/( (w - y)'*M*(w - y) + lambda * (v-u)'*M*(v-u));

	s = min(1, max(0, s));

	% Do the step
	u = u + s * (v - u);
	y = y + s * (w - y);

	% Plot solution
	U = zeros(n+1, n+1);
	U(2:end-1, 2:end-1) = reshape(u, n-1, n-1);
	surf(X1,X2,U);
	drawnow
end
