%Conferinta Brasov
p=1;
q=1;
a=0;b=1;
rT=0;
l1=-50
l=500;
Br=0.01; 

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


U=zeros(N,1);Uold=zeros(N+1,1);Uold(1)=1;
T=zeros(N,1);Told=zeros(N,1);
for i=1:N
    Told(i)=(1-rT)/(b-a)*(i-1)*h+(b*rT-a*1)/(b-a);
end

A1=zeros(N);
r1=zeros(N,1);
A2=zeros(N+1);
r2=zeros(N+1,1);


Start=1;
nr_it=0;
while Start==1
    
    %ecuatia temperaturii 
    %completam elementele matricei
    %bc1
    A1(1,1)=1;r1(1)=rT;
    %matricea
    for i=2:N-1
        A1(i,i-1)=1;A1(i,i)=-2;A1(i,i+1)=1;r1(i)=-Br*(Uold(i)-Uold(i-1))^2;
    end
    %bc2
    A1(N,N)=1;r1(N)=1;
    
    T=A1\r1;
    %ecuatia momentului 
    %completam elementele matricei
    %bc1
    A2(1,1)=1;r2(1)=0;
    %matricea
    for i=2:N-1
        A2(i,i-1)=1;A2(i,i)=-2;A2(i,i+1)=1;A2(i,N+1)=-h^2;r2(i)=-h^2*(l*T(i)+l1*T(i)^2);
    end
    %bc2
    A2(N,N)=1;r2(N)=0;
   
    %completam matricea cu ecuatia integralei
    for i=2:N-1
        A2(N+1,i)=1;
    end
    r2(N+1)=1/h;
    
    U=A2\r2;

    alfa=U(N+1);
    errU=norm(U(1:N)-Uold(1:N))/norm(U(1:N));
    errT=norm(T-Told)/norm(T);
  
    if (errU<1e-6)&(errT<1e-6)
        Start=2;
    else
            Uold=Uold+p*(U-Uold);
            Told=Told+q*(T-Told);
    end
    nr_it=nr_it+1;
    if mod(nr_it,1)==0
        disp([nr_it, errU, errT ])
    end
end
alfa
subplot(2,2,1)
plot(a:h:b,U(1:N),'b');
xlabel('Y');ylabel('U(Y)');title('Viteza')

subplot(2,2,2)
 plot(a:h:b,T,'r');
 xlabel('Y');ylabel('\theta(Y)');title('Temperatura')

subplot(2,2,3) 
y=a:h:b;
u=l.*(1-rT).*(-y.^3/6+y.^2/4-y/12)-6*y.^2+6*y;
plot(y,u,'k');
xlabel('Y');ylabel('U(Y)');title('Viteza pentru Br = 0 si {\lambda}^* = 0')

subplot(2,2,4)
plot(y,(1-rT)*y+rT,'k');
xlabel('Y');ylabel('\theta(Y)');title('Viteza pentru Br = 0 si {\lambda}^* = 0')




