function w = NonlinearFD()
                                                                                                 
alpha = 17;
beta = 43/3;
xmin = 1;
xmax = 3;
h = 0.1;
x = xmin:h:xmax;
M = 1000;
TOL = 1e-8;
k = 1;
                                                                                                 
n = length(x);
w = zeros(n,1);
v = zeros(n,1);
a = zeros(n,1);
b = zeros(n,1);
c = zeros(n,1);
d = zeros(n,1);
                                                                                                 
w(1) = alpha;
w(n) = beta;
                                                                                                 
for i = 2:n-1;
    w(i) = alpha+i*h*(beta-alpha)/(xmax-xmin);
end
                                                                                                 
v = w;
v(1) = 0;
v(n) = 0;
                                                                                                 
while (k <= M)
                                                                                                 
for i = 2:n-1;
    a(i) = 2+h^2*fy(x(i),w(i),(w(i+1)-w(i-1))/(2*h));
    b(i) =-1+h/2*fy1(x(i),w(i),(w(i+1)-w(i-1))/(2*h));
    c(i) = -1 -h/2*fy1(x(i),w(i),(w(i+1)-w(i-1))/(2*h));
    d(i) =-(2*w(i)-w(i+1)-w(i-1)+h^2*func(x(i),w(i),(w(i+1)-w(i-1))/(2*h)));
end
                                                                                                 
for i = 3:n-1;
    a(i) = a(i)-c(i)/a(i-1)*b(i-1);
    d(i) = d(i)-c(i)/a(i-1)*d(i-1);
end
                                                                                                 
v(n-1) = d(n-1)/a(n-1);
                                                                                                 
for i = n-2:-1:2;
    v(i) = (d(i) - v(i+1)*b(i))/a(i);
end
w = w + v;
                                                                                                 
if (norm(v,inf) < TOL)
    return;
end
                                                                                                 
k = k+1;
end
                                                                                                 
function y = func(x,y,y1)
y = 1/8*(32+2*x^3-y*y1);
                                                                                                 
function y = fy(x,y,y1)
y = -1/8*y1;
                                                                                                 
function y = fy1(x,y,y1)
y = -1/8*y;
                                                                                                 

