clear; close all; global epsilon; a = 1/sqrt(2); f1 = @(x) 1-exp(-(x(:,1)-a).^2-(x(:,2)-a).^2); global f2; f2 = @(x) 1-exp(-(x(:,1)+a).^2-(x(:,2)+a).^2); x1opts = []; x2opts = []; A = []; b = []; Aeq = []; beq = []; x0 = [0,0]; lb = [-4,-4]; ub = [4,4]; epsilons = linspace(0,1,20); for idx = 1:length(epsilons) epsilon = epsilons(idx); xopt = fmincon(f1, x0, A, b, Aeq, beq, lb, ub, @nonlcon); x1opts(end+1) = xopt(1); x2opts(end+1) = xopt(2); end xopts = [x1opts;x2opts]'; f1opts = f1(xopts); f2opts = f2(xopts); plot(f1opts,f2opts,'o'); hold on; x1 = linspace(0,1,20); y2 = interp1(f2opts, f1opts, x1,'spline'); hold on; plot(x1,y2); xlim([0,1]); ylim([0,1]); big; function [c,ceq] = nonlcon(x) global f2 epsilon; c = [f2(x)-epsilon]; ceq = []; end