s = linspace(-2, 2, 50);
e = ones(size(s));
lw = 'linewidth';

[x,y] = meshgrid(s, s);

f = @(x,y) exp(-(x.^2+y.^2));
fx = @(x,y) -2.*x.*f(x, y);
fy = @(x,y) -2.*y.*f(x, y);

z = f(x,y);

hf = figure(1); clf;

surf(x, y, z);
set(get(hf, 'children'), 'fontsize', 18);

shading flat
%camlight left

xlabel('x'); ylabel('y'); zlabel('z');
hold on;

axis equal

% choose a point
a = -0.25;
b = -0.1;
plot3(a, b, f(a,b), 'k*', lw, 4)
title(sprintf('a = %g, b = %g', a, b));

pause

plot3(a*e, s, f(a*e, s), 'r-', lw, 4);
plot3(s, b*e, f(s, b*e), 'r-', lw, 4);

camlight left

pause
s2 = linspace(a-1/2, a+1/2, 20);
e2 = ones(size(s2));
plot3(s2, b*e2, f(a,b) + fx(a,b).*(s2-a), 'b-', lw, 5);
plot3(a*e2, s2, f(a,b) + fy(a,b).*(s2-b), 'b-', lw, 5);
%plot3(s, b*e, fy(s, b*e), 'b-', lw, 2);


pause
[x2,y2] = meshgrid(linspace(a-1, a+1, 20), linspace(b-1, b+1, 20));

z2 = f(a,b) + fx(a,b).*(x2-a) + fy(a,b).*(y2-a);

mesh(x2, y2, z2)
