clf
[x,y]=meshgrid(linspace(1,3,200),linspace(1,3,201));
fun=@(x,y) (1+x).*(x-5).*(x-7).*sqrt(100-y.^4)/100;
mesh(x,y,fun(x,y)), pause
hold on
X=linspace(1.5,2.5,100);B=1.5*ones(size(X));
plot3([1.5 1.5],[1.5 1.5],[0,fun(1.5,1.5)],'y'), pause
plot(X,B,'k'), pause
plot3(X,B,fun(X,B),'k*'), pause
plot(B,X,'b'), pause
plot3(B,X,fun(B,X),'b*'), pause
f1=-.316667;
f2 = -.333392;

plot3(X,B,fun(1.5,1.5)+f1*(X-1.5),'k'), pause
plot3(B,X,fun(1.5,1.5)+f1*(X-1.5),'b'), pause
[x1,y1]=meshgrid(X,X);
z1=@(x,y) fun(1.5,1.5) + f1*(x1-1.5) + f2*(y1-1.5); 
mesh(x1,y1,z1(x1,y1)), pause
hold off
