Implementing a boundary value problem with code from file exchange

By the help of Ernesto Momox B. University of Essex, UK in Matlab file exchage, I got code for finite difference method that solves any boundary value problem. I am facing challeges editing the code to suite my taste. For example in the code, you have to create a separate m.file for your function and the input variables in the function handler is a bit different from what I would like it to be. I would like to edit the above code to solve the problem shown below with alpha = -1, beta = 1 and N=10 such that the code can take any real value of N. I would also like my function handle look likes this: nonlinearBVP_FDM(N, func, alpha, beta) with the fucntion:

But whenever I try to edit and run the code, matlab spews error messages in my command window as shown below:

% Index exceeds the number of array elements (10).
% Error in systemNxN (line 11)
% -w(9)+2*w(10)-w(11)+(h^2)*f(x(10),w(10),(w(11)-w(9))/(2*h));
% 
% Error in nonlinearBVP_FDM>@(w)systemNxN(w,x,h,alpha,beta) (line 50)
% w = fsolve(@(w) systemNxN(w,x,h,alpha,beta),w,options); % Solves the N x N
% 
% Error in fsolve (line 248)
%             fuser = feval(funfcn{3},x,varargin{:});
% 
% Error in nonlinearBVP_FDM (line 50)
% w = fsolve(@(w) systemNxN(w,x,h,alpha,beta),w,options); % Solves the N x N
% 
% Caused by:
%     Failure in initial objective function evaluation. FSOLVE cannot continue.

I post the code here for possoble addition to it and I will be glad. Thanks in advance.

function F = nonlinearBVP_FDM(a,b,alpha,beta)
% Nonlinear finite difference method for the general nonlinear 
% boundary-value problem -------------------------------------
% y''=f(x,y,y'), for a<x<b where y(a)=alpha and y(b)=beta. 
% ------------------------------------------------------------
% The interval [a,b] is divided into (N+1) equal subintervals
% with endpoints at x(i)=a+i*h for i=0,1,2,...,N+1. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Remarks: The function f should be defined as an m-file. %%
%%%          There is NO need for partial derivatives of f  %%
%%%          See given example                              %%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Example
% Solve the nonlinear boundary value problem
% y''=(1/8)*(32+2x^3-yy'), for 1<x<3, where y(1)=17 and y(3)=43/3
% Step 1... 
% Create the function f as a separate m-file and save it in the 
% current working directory.
% function f = f(x,y,yp)
% f = (1/8)*(32+2*x^3-y*yp); %Note that yp=y'
% Step 2... 
% In the command window type >> Y = nonlinearBVP_FDM(1,3,17,43/3);
% Note that Y(:,1) represents x and Y(:,2) is vector y(x) 
% The solution is then plotted in a new figure
% If the exact solution is given, plot it for comparison
% >> yexact = (Y(:,1)).^2+16./Y(:,1); plot(Y(:,1),yexact,'c')

format long
N=39;%Number of mesh points           
h=(b-a)/(N+1); %Step size
i=1:N;
x=a+i*h;
W=zeros(N+2,1); %Preallocate W and X
X=W;
w=alpha+i'*(beta-alpha)/(b-a)*h; %Initial approximations are obtained
% by passing a straight line through (a,alpha) and (b,beta)
options=optimset('TolFun',1E-10,'TolX',1E-10,'MaxIter',1E3,...
        'MaxFunEvals',1E3,'Display','off');

w = fsolve(@(w) systemNxN(w,x,h,alpha,beta),w,options); % Solves the N x N
% nonlinear system of equations

W(1)=alpha;         X(1)=a;
W(2:length(w)+1)=w; X(2:length(x)+1)=x;
W(N+2)=beta;        X(N+2)=b;

F = [X W];

hand = figure;
set(hand,'Name','Nonlinear Finite Difference Method','NumberTitle','off');
plot(X,W,'r.','MarkerSize',16), grid on, hold on
plot(X(1),W(1),'bo','MarkerSize',10,'LineWidth',2)
plot(X(N+2),W(N+2),'go','MarkerSize',10,'LineWidth',2)
xlabel('$x$','FontSize',20,'interpreter','latex')
ylabel('$y(x)$','FontSize',20,'interpreter','latex')
legend('Numerical Solution','y(a) = alpha','y(b) = beta','Orientation',...
       'horizontal','Location','Best')