% function y = nonlinshoot(a,b,n,alpha,beta,f,fy,fyp,tol,maxit)
%
% Non-linear finite differences to solve BVP
%
%   y'' = f(t,y,y') 
%   y(a) = alpha
%   y(b) = beta
%
%  using Newton's method and variational equation
% 
% Inputs
%  a,b    time interval
%
%  n      number of time subintervals
%
%  alpha  = y(a)
%
%  beta   = y(b)
%
%  f      function handle, for a general non-linear problem
%         example:
%
%         f = @(t,y,yp) -3*yp +2*y+2*t+3; 
%
%         the parameters are (in this order)
%         t   space variable
%         y   y(t)
%         yp  y'(t)
%
%  fy     function handle for derivative of f with respect to
%         y argument (SECOND argument)
%
%  fyp    function handle for derivative of f with respect to
%         yp argument (THIRD argument)
%
%  tol    tolerance for stopping iterations
%
%  maxit  maximum number of iterations
%
% Outputs
%  y      function y at the points t = a + i*(b-a)/n, i=0..n.

function y = nonlinfd(a,b,n,alpha,beta,f,fy,fyp,tol,maxit)
% setup grid
h=(b-a)/(n+1);  t = a + (0:n+1)'*h;

% initial guess is line satisfying bc
y = % TO BE COMPLETED

% discrete Laplacian (for later)
L = (1/h^2)*spdiags( TO BE COMPLETED )


for k=1:maxit,
 % evaluate non-linear map F

 % output
 fprintf('it=%d, norm(F)=%e\n',nit,norm(F));

 % check convergence
 if (norm(F)<tol)
  return;
 end;

 % build Jacobian of F
 DF = sparse(n,n);
 To BE COMPLETED

 % Newton update (we update only the parts of y that are not BC)
 y(2:end-1) = y(2:end-1) + TO BE COMPLETED
end;

fprintf('max number of iterations (%d) reached and error is %g\n',maxit,norm(F));
