> restart:with(VectorCalculus):with(plots):with(plottools):z:=(x,y)->sin(1/2*x^2-1/4*y^2+3)*cos(2*x+1-exp(y)):
> z:=(x,y)->x+y/2:
> t1:=(u,v,w)->(exp(u+v),exp(u-v),cos(w));
> ti:=[t1(x,y,z)];
> bt:=[x,y,z(x,y)];
> bl:=t1(bt[]);
> pc:=[1/2,1/2];
> p0:=[pc[1],pc[2],z(pc[1],pc[2])];
> myunit:=0.3;
> circum:=map(p->[p[1],p[2],z(p[1],p[2])],[(pc+[myunit,i/5*myunit])$i=-10..10,(pc+[-myunit,i/10*myunit])$i=-10..10,(pc+[i/10*myunit,myunit])$i=-10..10,(pc+[i/10*myunit,-myunit])$i=-10..10]):
> i0:=evalf([t1(p0[])]);
> display(plot3d([x,y,z(x,y)],x=-1.2..1.2,y=-1.2..1.2,axes=normal),point(p0,color=blue,symbol=circle,symbolsize=10),'point(circum[i],color=red,symbol=circle,symbolsize=2)'$'i'=1..nops(circum));
> display(plot3d([bl],x=-1.2..1.2,y=-1.2..1.2,axes=normal),point(i0,color=blue,symbol=circle,symbolsize=10),'point([t1(circum[i][])],color=red,symbol=circle,symbolsize=2)'$'i'=1..nops(circum));
> J:=Jacobian(ti,[x,y,z]);
> lindiff:=Matrix(eval(J,{x=p0[1],y=p0[2],z=p0[3]}));
> tl:=(u,v,w)->convert(lindiff.Vector([u,v,w]),list);
> evalf(tl(0,0,0));
> evalf([t1(p0[])]);
> linapp:=map(c->evalf(i0+tl((c-p0)[])),circum):
> display(plot3d([bl],x=-1.2..1.2,y=-1.2..1.2,axes=normal),point(i0,color=blue,symbol=circle,symbolsize=10),'point([t1(circum[i][])],color=red,symbol=circle,symbolsize=2)'$'i'=3..nops(circum),'point(linapp[i],color=green,symbol=circle,symbolsize=2)'$'i'=4..nops(circum));
> display('point(linapp[i]-[t1(circum[i][])],color=green,symbol=circle,symbolsize=2)'$'i'=1..nops(circum));