dx = 0.01;
z01 = 5;
z02 = 100;
X1 = [T0;z01];
diff_eq1 = @(x,X1) [X1(2), a*((Ta+Tb)^4 + 4*(Ta+Tb)^3*(X1(1)-Tb)) + b*(X1(1)-Tb)];
[~,Y1] = rk4(diff_eq1,dx,X1,L,x0);
function [time,y] = rk4(diff_eq,h,y0,T,t0)
y(:,1) = y0;
t(1) = t0;
tf = T;
time=t(1):h:tf;
for i = 1:length(time)-1
k1(:,i) = diff_eq(t(i),y(:,i));
k2(:,i) = diff_eq(t(i)+h/2,y(:,i)+k1*h/2);
k3(:,i) = diff_eq(t(i)+h/2,y(:,i)+k2*h/2);
k4(:,i) = diff_eq(t(i)+h,y(:,i)+k3*h);
y(:,i+1) = y(:,i) + (1/6)*(k1(:,i) + 2*k2(:,i) + 2*k3(:,i) + k4(:,i))*h;
t(i+1) = t(i)+h;
y;
end
end