function bvpscripts
h=pi/2/4;
A=[1 0 0 0 0; 1 -(2-h^2) 1 0 0; 0 1 -(2-h^2) 1 0; ...
    0 0 1 -(2-h^2) 1; 0 0 0 0 1];
b=[1; h^2; h^2; h^2; 0];
x=linspace(0,pi/2,5);
y=A\b;
plot(x,y)
disp('press any key to continue')
pause

N=20;
h=pi/2/N;
A=-(2-h^2)*eye(N+1);
A(1,1)=1;
A(N+1,N+1)=1;
for i=1:N-1
    A(i+1,i)=1;
    A(i+1,i+2)=1;
end
b=h^2*ones(N+1,1);
b(1)=1;
b(N+1)=0;
x=linspace(0,pi/2,N+1);
y=A\b;
plot(x,y)
disp('press any key to continue')
pause

N=20;
h=pi/2/N;
v=ones(N,1);
A=-(2-h^2)*eye(N+1)+diag(v,1)+diag(v,-1);
A(1,2)=0; A(N+1,N)=0;
A(1,1)=1; A(N+1,N+1)=1;
b=h^2*ones(N+1,1);
b(1)=1;
b(N+1)=0;
x=linspace(0,pi/2,N+1);
y=A\b;
plot(x,y)
xlabel('x')
ylabel('y')
disp('press any key to continue')
pause

xlow=0;
xhigh=pi/2;
solinit = bvpinit(linspace(xlow,xhigh,10),[1 -1]);
sol = bvp4c(@bvp4ode,@bvp4bc,solinit);
xint = linspace(xlow,xhigh);
Sxint = deval(sol,xint);
plot(xint,Sxint(1,:))
xlabel('x')
ylabel('y')
disp('press any key to continue')
pause

h=1/4
A=[1 0 0 0 0; 
   1 -(2-h^2) 1 0 0; 
   0 1 -(2-h^2) 1 0;
   0 0 1 -(2-h^2) 1;
   0 0 0 2 -(2-h^2)];
B=[1; 0; 0; 0; 0];
y=A\B;
plot(0:h:1,y)
disp('press any key to continue')
pause

xlow=0;
xhigh=1;
solinit = bvpinit(linspace(xlow,xhigh,10),[1 -1]);
sol = bvp4c(@bvp5ode,@bvp5bc,solinit);
xint = linspace(xlow,xhigh);
Sxint = deval(sol,xint);
plot(xint,Sxint(1,:))
xlabel('x')
ylabel('y')

% -----------------------------------------------
function dydx = bvp4ode(x,y,lambda)
dydx = [  y(2)
         1-y(1) ];
% -----------------------------------------------
function res = bvp4bc(ya,yb,lambda)
res = [  ya(1)-1
         yb(1) ];
%
function dydx = bvp5ode(x,y)
dydx = [  y(2)
         -y(1) ];
% -----------------------------------------------
function res = bvp5bc(ya,yb)
res = [  ya(1)-1
         yb(2) ];