

a=0;b=1;
L=500;%lambda
rT=0.5;
G=1;%gamma
D=1;%Darcy 

N=101;
h=(b-a)/(N-1);

U=zeros(N,1);
T=zeros(N,1);
for i=1:N
    y=(i-1)*h;
    T(i)=-0.5*y^2+(2*rT+0.5)*y-rT;
end
figure(1)
plot(a:h:b,T)


A=zeros(N+1);
r=zeros(N+1,1);


 
    %ecuatia momentului 
    %completam elementele matricei
    %bc1
    A(1,1)=1;r(1)=0;
    %matricea
    for i=2:N-1
        A(i,i-1)=1;A(i,i)=-(2+h*h/G/D);A(i,i+1)=1;A(i,N+1)=-h^2/G;r(i)=-L*h^2/G*T(i);
    end
    %bc2
    A(N,N)=1;r(N)=0;
   
    %completam matricea cu ecuatia integralei
    for i=2:N-1
        A(N+1,i)=1;
    end
    r(N+1)=1/h;
    
    U=A\r;

    alfa=U(N+1);


alfa
figure(2)
plot(a:h:b,U(1:N),'b');
