clear, clc
L = 10;
solinit = bvpinit(linspace(0,L),[1200 0]); 
options = bvpset('RelTol', 1e-6);
sol = bvp4c('RK','bc2',solinit,options); 
plot(sol.x, sol.y(1,:)), 
grid on 
legend('T(x)','dT(x)') 
xlabel('x')

clear, clc
L = 10;
solinit = bvpinit(linspace(0,L),[1200 0]); 
options = bvpset('RelTol', 1e-6);
sol = bvp4c('RK','bc2',solinit,options); 
plot(sol.x, sol.y(1,:)), 
grid on 
legend('T(x)','dT(x)') 
xlabel('x')

function dTETA=RK(x,TETA)
dTETA=zeros(2,1);
Y = TETA(2); % производная
alw=50; 
Rw=0.05; 
R1=0.3; 
R2=1500;
c2=660; 
la2=4; 
la1=0.084; 
al1=0.5; 
al2=0.5;
et0=0.077; 
R=8.314;
cp=1000; 
Tv=315;
T0=300; 
E=12; 
Q=1513000;
k0=50000000000;
u=0.001;
v=0.5;
u0=v/u;
L=0.5;
Te=T0+(Q*et0)/(1+(R2*c2*u)/(R1*cp*(v+u)));
%TETA=(Tv-T0)/(Te-T0);
Al0=(2*alw)/Rw;
a1=(R1*(u0+1)*cp+R2*c2*u0)/(al1*la1+al2*la2);
a2=Al0/(al1*la1+al2*la2);
a3=(R1*Q*et0*k0)/(al1*la1+al2*la2);
dTETA(1)=Y; % первая производная dT/dx
dTETA(2)=L*u*a1*Y+L*L*a2*TETA(1)-((a3*L*L)/(T0+TETA(1)*(Te-T0)))*exp(-E/(R*(T0+TETA(1)*(Te-T0))));
end

function res=bc2(ua,ub) % граничные условия 
% T -> u_(1)
% dT -> u_(2)
% ua -- x==0
% ub -- x==L 
res = [ua(1)-1200; % T(0)=1200
   ub(1)-300]; % T(L)=300
end

function dTETA=RK2(x,TETA)
dTETA=zeros(2,1);
 
dTETA(1) = TETA(2); % первая производная dT/dx
dTETA(2) = TETA(2) + TETA(1) - 1/(TETA(1))*exp(-1/TETA(1));
end

function res=bc2(ua,ub) % граничные условия 
res = [ua(1)-1200; % T(0)=1200
       ub(1)-300]; % T(L)=300
end

clear all
clc
L = 10;
solinit = bvpinit([0 L],[1200 0]); 
options = bvpset('RelTol', 1e-6,'NMax',5e4);
sol = bvp4c(@RK2,@bc2,solinit,options); 
plot(sol.x, sol.y(1,:)), 
grid on 
%legend('T(x)') 
xlabel('x')

function [  ] = Temp( T1,T2,T0,A,L )
 
s=figure;
step=0.1;
t=0;
h=1;
Ti=T0*ones(L,1);
Tk=T0*ones(L,1);
colormap(bone);
 
while ishandle(s)
    Ti(1)=T1;
    Ti(L)=T2;
    Tk(1)=T1;
    Tk(L)=T2;
    for x=2:h:(L-1)
        Tk(x)=(A*step/h^2)*(Ti(x-1)+Ti(x+1)-2*Ti(x)) + Ti(x);
    end
    image(Ti);
    Ti=Tk;
    t=t+step;
    pause(step);
end
 
end

>> N =10; % количество узлов в сетке;
t_end = 60; %окончание по времени;
L = 0.3; %длина толщены пластинки;
T0 = 300; %начальная температура;
T1 = 300; % температура на границе x=0);
Tr = 100; %температура на границе x=L ;
alw=50; 
Rw=0.05; 
R1=0.3; 
R2=1500;
c2=660; 
la2=4; 
la1=0.084; 
al1=0.5; 
al2=0.5;
et0=0.077; 
R=8.314;
cp=1000; 
Tv=315;
E=126000; 
Q=1513000;
k0=50000000000;
u=0.001;
v=0.5;
u0=v/u;
Al0=(2*alw)/Rw;
Fi=(R2*c2)/(R1*cp);
kap=(al1*la1+al2*la2)/(R2*c2+R1*cp);
a0=Al0/(R2*c2+R1*cp);
a1=(R1*T0*k0*Q*et0)/(R2*c2+R1*cp);
a2=v/(1+Fi);
T = zeros(1,N);
alfa = zeros(1,N);
beta = zeros(1,N);
 
h = L/(N-1); % определяем расчетный шаг сетки по пространственной координате
tau = t_end/100; %определяем расчетный шаг сетки по времени
for i=1:N
    T(i)= T0;
    time = 0;
    
    while time < t_end
        time = time + tau;
        alfa(1) = 0;
        beta(1) = 1;
        
        for j=2:N-1
            ai = kap/(h^2);
            bi = 1/tau+(2*kap)/(h^2)+a2/h;
            ci = a2/h+kap/(h^2);
            fi = -1/tau+a0*(T(i)-T0)-a1/(T(i))*exp(-E/R/T(i));
            alfa(j) = ai/(bi-ci*alfa(j-1)); % прогоночные коэффициенты
            beta(j) = (ci*beta(j-1)-fi)/(bi-ci*alfa(j-1));% прогоночные коэффициенты
        end
        
        T(N) = Tr; %определяем значение температуры на правой границе
        
        for k = N-1:-1:1
            T(k)= alfa(k)*T(k+1)+ beta(k);
        end
    end
end
plot(h*(0:N-1),T);