clc;
clear;
global h X  d
a=0;
b=1;
m=10; 
N=100;
%k=(b-a)/N;
h=(b-a)/m;

nn=100;
hh=3*h/nn;
X=-2*h:h:1+2*h;
X=X';

U=(X(3:m+3)).^2.*(1-(X(3:m+3))).^2;
B=zeros(m-2,m-2);
C=zeros(m-2,1);
 d=zeros(m+2,N+2);

 

for j=1:m-2;
    for i=1:m-2;
        fun1=@(x) (((i+1).*h-x).^2-3*((i).*h-x).^2+3*((i-1).*h-x).^2).*(((i+1).*h-x).^2-3*((i).*h-x).^2+3*((i-1)*h-x).^2);
        fun2=@(x) (((i+1).*h-x).^2-3*((i).*h-x).^2).*(((i+1).*h-x).^2-3*((i).*h-x).^2);
        fun3=@(x) (((i+1).*h-x).^2).*(((i+1).*h-x).^2);
        fun4=@(x) (((i+1).*h-x).^2).*(((j+1).*h-x).^2-3*((j).*h-x).^2);
        fun5=@(x) (((i+1).*h-x).^2-3*((i).*h-x).^2).*(((j+1).*h-x).^2-3*((j).*h-x).^2+3*((j-1).*h-x).^2);
        fun6=@(x) (((i+1).*h-x).^2).*(((j+1).*h-x).^2-3*((j).*h-x).^2+3*((j-1).*h-x).^2);
        fun7=@(x) (((j+1).*h-x).^2).*(((i+1).*h-x).^2-3*((i).*h-x).^2);
        fun8=@(x) (((j+1).*h-x).^2-3*((j).*h-x).^2).*(((i+1).*h-x).^2-3*((i).*h-x).^2+3*((i-1).*h-x).^2);
        fun9=@(x) (((j+1).*h-x).^2).*(((i+1).*h-x).^2-3*((i)*h-x).^2+3*((i-1).*h-x).^2);
        fun10=@(x) (((i+1).*h-x).^2).*(x.^2.*(1-x).^2);
        fun11=@(x) (((i+1).*h-x).^2-3*((i).*h-x).^2).*(x.^2.*(1-x).^2);
        fun12=@(x) (((i+1).*h-x).^2-3*((i).*h-x).^2+3*((i-1).*h-x).^2).*(x.^2.*(1-x).^2);
        if (i==j)
            if(i==1) 
               B(1,1)=(1/h^4)*integral(fun3,0,h,'RelTol',1e-12,'AbsTol',1e-13);
            elseif(i==2)
               B(i,j)=(1/h^4)*(integral(fun2,0,h,'RelTol',1e-12,'AbsTol',1e-13)+integral(fun3,h,2*h,'RelTol',1e-12,'AbsTol',1e-13));
            elseif(i>=3 && i<=m)
               B(i,j)=(1/h^4)*(integral(fun1,(i-3)*h,(i-2)*h,'RelTol',1e-12,'AbsTol',1e-13)+integral(fun2,(i-2)*h,(i-1)*h,'RelTol',1e-12,'AbsTol',1e-13)+integral(fun3,(i-1)*h,i*h,'RelTol',1e-12,'AbsTol',1e-13));
            elseif(i==m+1)
               B(i,j)=(1/h^4)*(integral(fun1,(i-3)*h,(i-2)*h,'RelTol',1e-12,'AbsTol',1e-13)+integral(fun2,(i-2)*h,(i-1)*h,'RelTol',1e-12,'AbsTol',1e-13));
            else
               B(i,j)=(1/h^4)*(integral(fun1,(i-3)*h,(i-2)*h,'RelTol',1e-12,'AbsTol',1e-13));
            end
        end
        
        if (i==j+1)
            if(j==1)
             B(i,j)=(1/h^4)*integral(fun7,0,h,'RelTol',1e-12,'AbsTol',1e-13);
            elseif(i>=3 && i<=m+1)
                B(i,j)=(1/h^4)*(integral(fun8,(j-2)*h,(j-1)*h,'RelTol',1e-12,'AbsTol',1e-13)+integral(fun7,(j-1)*h,j*h,'RelTol',1e-12,'AbsTol',1e-13));
            else
                B(i,j)=(1/h^4)*integral(fun8,(j-2)*h,(j-1)*h,'RelTol',1e-12,'AbsTol',1e-13);
            end
        end
        if (i==j+2)
            B(i,j)=(1/h^4)*integral(fun9,(j-1)*h,j*h,'RelTol',1e-12,'AbsTol',1e-13);
        end
      
        if (j==i+1)
           if(i==1)
            B(i,j)=(1/h^4)*(integral(fun4,0,h,'RelTol',1e-12,'AbsTol',1e-13));
           elseif(i>=2 && i<=m)
             B(i,j)=(1/h^4)*(integral(fun5,(i-2)*h,(i-1)*h,'RelTol',1e-12,'AbsTol',1e-13)+integral(fun4,(i-1)*h,i*h,'RelTol',1e-12,'AbsTol',1e-13));
           else
             B(i,j)=(1/h^4)*(integral(fun5,(i-2)*h,(i-1)*h,'RelTol',1e-12,'AbsTol',1e-13));
            end
        end
        if (j==i+2)
            B(i,j)=(1/h^4)*(integral(fun6,(i-1)*h,i*h,'RelTol',1e-12,'AbsTol',1e-13));
        end
       if(i==1)
            C(i)=(1/h^2)*integral(fun10,0,h,'RelTol',1e-12,'AbsTol',1e-13);
        elseif(i==2)
            C(i)=(1/h^2)*(integral(fun11,0,h)+integral(fun10,h,2*h,'RelTol',1e-12,'AbsTol',1e-13));
        elseif(i>=3 && i<=m)
            C(i)=(1/h^2)*(integral(fun12,(i-3)*h,(i-2)*h,'RelTol',1e-12,'AbsTol',1e-13)+integral(fun11,(i-2)*h,(i-1)*h,'RelTol',1e-12,'AbsTol',1e-13)+integral(fun10,(i-1)*h,i*h,'RelTol',1e-12,'AbsTol',1e-13));
        elseif(i==m+1)
            C(i)=(1/h^2)*(integral(fun12,(i-3)*h,(i-2)*h,'RelTol',1e-12,'AbsTol',1e-13)+integral(fun11,(i-2)*h,(i-1)*h,'RelTol',1e-12,'AbsTol',1e-13));
        else
            C(i)=(1/h^2)*integral(fun12,(i-3)*h,(i-2)*h,'RelTol',1e-12,'AbsTol',1e-13);
            
        end
    end
end
d(3:m,1)=B\C;
F=zeros(m+1,1);
S=zeros(m+1,1);

for j=3:m+3
    for i=1:m+2
        F(j-2)=F(j-2)+d(i,1)*bsp2(X(j),i);
        
% % z=X(j):hh:X(j+2);
% % s=size(z,2);
% % for t=1:s
% %     y(t)=bsp2(z(t),i);
% % end
% % plot(z,y)
% % hold on
% % for ii=1:10000
% %     disp('')
% % end
    end
    S(j-2)=abs((F(j-2)-U(j-2)));
end

 p=norm(S,Inf);
    plot(X(3:m+3),F);
    hold on
    plot(X(3:m+3),U);
 p