setup;

% Initial guess
u = ua;

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

% Line search parameter
gamma = 0.1;

for i=1:100
	% Solve adjoint equation K*p = y - yd
	p = K \ (y - yd);

	% Evaluate objective
	obj = 1/2*(y - yd)'*M*(y-yd) + lambda/2*u'*M*u;

	% Compute search direction
	v = -(beta.*p + lambda*u);

	% Compute step length s
	s = 1/lambda;
	while true

		u_step = min(ub, max(ua, u + v * s));
		y_step = K \ (beta .* u_step);

		obj_step = 1/2*(y_step - yd)'*M*(y_step-yd) + lambda/2*u_step'*M*u_step;

		if obj_step <= obj - gamma/s *(u_step - u)'*M*(u_step-u);
			break
		end

		s = s / 2;

		if s < 1e-10
			error('No valid step length');
		end

	end

	s

	% Do the step
	u = u_step;
	y = y_step;

	% 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
