nlbvpfd2

% Program 7.1 Nonlinear Finite Difference Method for BVP
% Uses Multivariate Newton?s Method to solve nonlinear equation
% Inputs: interval inter, boundary values bv, number of steps n
% Output: solution w
% Example usage: W=nlbvpfd([0 1], [1 4], 40)
function w=nlbvpfd2(inter,bv,n)
x = linspace(0, 1, n);
a=inter(1); b=inter(2); ya=bv(1); yb=bv(2);
h=(b-a)/(n+1); % h is step size
w=zeros(n,1); % initialize solution array w
for i=1:20 % loop of Newton step
  w=w-jac(w,inter,bv,n)\f(w,inter,bv,n);
end
plot([a a+(1:n)*h b],[ya w' yb]); % plot w with boundary data

function y=f(w,inter,bv,n)
x = linspace(0, 1, n);
y=zeros(n,1);h=(inter(2)-inter(1))/(n+1);
y(1)=w(2)-2*w(1)+x(1)^2*w(2)-2*x(1)^2*w(1)+x(1)^2-2*h*w(2)-2*h-6*w(1);
y(n)=(-4/3)-2*w(n)+w(n-1)+x(n)^2*(-4/3)-2*x(n)^2*w(n)+x(n)^2*w(n-1)-2*h*(-4/3)-2*h*w(n-1)-6*w(n);
for i=2:n-1
  y(i)=w(i+1)-2*w(i)+w(i-1)+x(i)^2*w(i+1)-2*x(i)^2*w(i)+x(i)^2*w(i-1)-2*h*w(i+1)-2*h*w(i-1)-6*w(i);
end


function a=jac(w,inter,bv,n)
x = linspace(0, 1, n);
a=zeros(n,n);h=(inter(2)-inter(1))/(n+1);
a(n,n) = -2-2*x(n)^2-6;
a(1,1) = -2-2*x(1)^2-6;
for i=2:n-1
  a(i,i)=-2-2*x(i)^2-6;
end
a(1,2) = 1+x(1)^2-2*h;
a(2,1) = 1+x(2)^2-2*h;
a(n-1,n) = 1+x(n-1)^2-2*h;
a(n,n-1) = 1+x(n)^2-2*h;
for i=2:n-2
  a(i,i+1)=1+x(i)^2-2*h;
  a(i+1,i)=1+x(i)^2-2*h;
end