function Y = blbvpfd(inter,bv,n)

a = inter(1);
b = inter(2);
y0 = bv(1);
yN = bv(2);
h = (b-a)/(n+1);
Y = zeros(n,1);
for i = 1:20
    Y = Y - DF(Y,inter,bv,n)\F(Y,inter,bv,n);
end
plot([a a+(1:n)*h b],[y0 Y' yN]);
end
function Ynew = F(Y,inter,bv,n)
Ynew = zeros(n,1);
h = (inter(2)-inter(1))/(n+1);
Ynew(1) = bv(1)-(2+h^2)*Y(1) + h^2*Y(1)^2+Y(2);
Ynew(n) = Y(n-1) -(2+h^2)*Y(n) + h^2*Y(n)^2+bv(2);
for i = 2:n-1
    Ynew(i) = Y(i-1) - (2+h^2)*Y(i)+h^2*Y(i)^2 + Y(i+1);
end

end
function a = DF(Y,inter,bv,n)
a = zeros(n,n);
h = (inter(2)-inter(1))/(n+1);
for i =1:n
    a(i,i) = 2*h^2*Y(i)-2-h^2;
end

for i =1:n-1
    a(i,i+1) = 1;
    a(i+1,i) = 1;
end
end
