# newton_sys.R -- R.Taylor 19.11.2018
# R code implementing Newton's method to solve the non-linear system
#   3 x^2 - y^2 = 0
#   2 x y^2 - x^3 - 1 = 0

# initial guess for solution:
# (this syntax constructs a vector of 2 components)
x <- c(0.8,1.0);

for (n in 1:5) {  # loop to do 5 iterations of Newton's method

  # evaluate the vector y=f(x)
  # (this sytax constructs a vector with two components)
  f <-  c( 3*x[1]^2 - x[2]^2, 2*x[1]*x[2]^2 - x[1]^3 - 1 );

  # evaluate matrix of partial derivatives Df(x):
  # (this syntax usies elements of a vector to fill the matrix column-by-column):
  D <- matrix( c( 6*x[1],
                  2*x[2]^2-3*x[1]^2,
                 -2*x[2],
                  4*x[1]*x[2] ),
               nrow=2, ncol=2 );

  # solve the linear system (D)(dx) = f for the vector dx:
  dx <- solve( D, f );

  # print some info about the current iteration:
  cat("\n== Iteration", n, "=======================================\n\n");
  cat("x =", x, "\n\n");
  cat("f(x) =",f, "\n\n" );
  cat("matrix of partial derivatives Df(x):\n")
  print(D);

  cat("\ndx =", dx,"\n\n");

  # update the solution using our calculated dx:
  x <- x - dx;

  cat("after update: x =", x, "\n");

}