N=10;
a=0;
b=0;
c=10000;
A=zeros(N);
B=zeros(N,1);

s=zeros(1,N+2);
s(1,N+2)=1;
s(1,2:N+1)=sort(rand(1,N));
y=zeros(1,N+2);

for i=2:N+1
    A(i-1,i-1)=c/3*(s(1,i+1)-s(1,i-1))-1/(s(1,i)-s(1,i-1))-1/(s(1,i+1)-s(1,i));
end

for i=2:N
    A(i,i-1)=1/6*c*(s(1,i+1)-s(1,i))+1/(s(1,i+1)-s(1,i));
    A(i-1,i)=A(i,i-1);
end

for i=2:N+1
    
    l=s(1,i-1);
    m=s(1,i);
    n=s(1,i+1);
    
    B(i-1,1)=m-l+1/3*c*(b-a)/(m-l)*(m^3-l^3)-1/2*c*(b-a)*l/(m-l)*(m^2-l^2)+1/2*c*a*(m-l)...
            +n-m+1/3*c*(b-a)/(n-m)*(n^3-m^3)-1/2*c*(b-a)*l/(n-m)*(n^2-m^2)+1/2*c*a*(n-m);
 end
 
 v=A\B;
 
v=v';
y(1,2:N+1)=v;
y=y+(b-a)*s+a;

plot(s,y);