d=[];
n=20;
x0=0; xend=2;
y0=0; yend=1;
h=(xend-x0)/n;
for  i=1:n-1
	d=[d,-2+i*h^3];
end;
l=ones(1,n-2);
u=l; 
b=zeros(size(d));
b(1)=-y0;
b(end)=-yend; 
y=trisys(u,d,l,b);

x=linspace(x0,xend,n+1);
y=[y0,y,yend]

plot(x,y,'o',x,x.*hyperg_0F1(4/3,-x.^3/9)/hyperg_0F1(4/3,-8/9)/2);

I=[1:n-1,2:n-1,1:n-2];
J=[1:n-1,1:n-2,2:n-1];
V=[-2+(1:n-1)*h^3,ones(1,2*n-4)];
A=sparse(I,J,V);
B=zeros(n-1,1);
B(1)=-y0;
B(end)=-yend;
Y=A\B;

Y=[y0,Y',yend];
y-Y
