% tangent plane approximation
% to a surface
close all

f = @(x,y) exp(1./(x.^2 + y.^2 + 1));
[xx,yy] = meshgrid(-2:0.05:2);
[xxp,yyp] = meshgrid(-1:0.05:1);
[fx,fy] = gradient(f(xx,yy),0.05);

x0 = 0.0;
y0 = 0.0;
t = (xx == x0) & (yy == y0);
indt = find(t);
fx0 = fx(indt);
fy0 = fy(indt);

z = @(x,y) f(x0,y0) + fx0*(x-x0) + fy0*(y-y0);

surf(xx,yy,f(xx,yy),'EdgeAlpha',0.7,'FaceAlpha',0.9)
hold on
surf(xxp,yyp,z(xxp,yyp))
plot3(x0,y0,f(x0,y0),'r*')