function TaylorExpasion
x0 = [-2, 2]
i = 0;
iMax = 1000;
x = x0;
Delta = 1;
alpha = 0.5;

figure; hold on
pnts = [];
while i < iMax && Delta>10^(-3)
    p = -hessian(x)^(-1).*jaco(x);       
    p = diag(p)';  
    xold = x;
    x = x + alpha.*p;
    Delta = sum( (x-xold).^2 );
    fprintf( '[%f, %f]\n', x(1), x(2) );
    pnts = [pnts; x];
    i = i+1;
end

[X,Y] = meshgrid(-2:.1:2, -2:.1:2); 
surf(X,Y,func(X,Y));
plot3( pnts(:,1), pnts(:,2), func(pnts(:,1), pnts(:,2)), '-sr' );
axis([-2 2 -2 2 0 2])

function v = func(x, y) % the function
v = 0.5*x.^2 + 0.2*y.^2;

function v = jaco(x) % the jacobian matrix of the first derivatives
v = [x(1), 0; 0, 0.4*x(2)];

function v = hessian(x) % the hessian matrix of second derivetives 
v = [1, 0; 0, 0.4];