The classic Newton method approximates the function to be optimised 
as a quadratic using the Taylor series expansion:
=f %2B \nabla f^T p %2B 1/2 p^T \nabla ^2fp)
By minimising this function with respect to 
the optimal search direction can be found as:
The step length 
is then computed via a backtrack linesearch using Wolfe conditions that assure sufficient decrease.
The inverse of the Hessian 
is computationally expensive to compute due to both finite difference limitations and the cost of inverting a particularly large matrix. For this reason an approximation to the inverse of the Hessian is used 
.
This approximation is updated at each iteration based on the change in 
and the change in 
as follows:
H(I-\frac{ys^T}{y^Ts}) %2B \frac{ss^T}{y^Ts})
For more details on implementation I highly advise Nocedal's book, Numerical Optimization.
http://www.apmath.spbu.ru/cnsa/pdf/monograf/Numerical_Optimization2006.pdf
Testing the BFGS algorithm on the Rosenbrock function in 2 dimensions, an optimal solution is found in 34 iterations.
The code implements an initial Hessian 
as the identity matrix, and if the problem is two dimensional then the code can produce a trajectory plot of the optimisation scheme.
The central difference method is used for the calculation of gradients.
