function thermal_boundary1
clc
clear all
close all
x=0;
x=4;
xl=fzero(@eqsol,x)
%%Parameters of Blasius Equation
U_inf = 1;
L = 10;
mu = 1.789E-5;
rho = 1.225;
nu = mu/rho;
A = sqrt(nu/U_inf);
h = 0.01;
%%Numerical Solution of Blasius Equation
dq1 = @(x, r1, r2, r3) r2;
dq2 = @(x, r1, r2, r3) r3;
dq3 = @(x, r1, r2, r3) -r1*r3;
eta = 0:h:10;
x = 0:h:10;
r1(1) = 0;
r2(1) = 0;
r3(1) = 0.4696;
for i = 1:(length(eta)-1)
a = h.*[dq1(eta(i), r1(i), r2(i), r3(i)), dq2(eta(i), r1(i), r2(i), r3(i)), dq3(eta(i), r1(i), r2(i), r3(i))];
b = h.*[dq1(eta(i), r1(i)+a(1)/2, r2(i)+a(2)/2, r3(i)+a(3)/2), dq2(eta(i)+h/2, r1(i)+a(1)/2, r2(i)+a(2)/2, r3(i)+a(3)/2), dq3(eta(i)+h/2, r1(i)+a(1)/2, r2(i)+a(2)/2, r3(i)+a(3)/2)];
c = h.*[dq1(eta(i), r1(i)+b(1)/2, r2(i)+b(2)/2, r3(i)+b(3)/2), dq2(eta(i)+h/2, r1(i)+b(1)/2, r2(i)+b(2)/2, r3(i)+b(3)/2), dq3(eta(i)+h/2, r1(i)+b(1)/2, r2(i)+b(2)/2, r3(i)+b(3)/2)];
d = h.*[dq1(eta(i), r1(i)+c(1), r2(i)+c(2), r3(i)+c(3)), dq2(eta(i)+h, r1(i)+c(1), r2(i)+c(2), r3(i)+c(3)), dq3(eta(i)+h, r1(i)+c(1), r2(i)+c(2), r3(i)+c(3))];
r3(i+1) = r3(i)+ 1/6*(a(3)+2*b(3)+2*c(3)+d(3));
r2(i+1) = r2(i)+ 1/6*(a(2)+2*b(2)+2*c(2)+d(2));
r1(i+1) = r1(i)+ 1/6*(a(1)+2*b(1)+2*c(1)+d(1));
end
%%Plotting and Visualization
figure(2)
plot(eta,r1,eta, r2, eta, r3, 'LineWidth', 2)
xlim([0 10])
title('Solution of Blasius eqution', 'FontSize', 14);
xlabel('f, f'' and f''''', 'FontSize', 20);
ylabel('\eta', 'FontSize', 20);
grid on
Legend1 = {'f(\eta)', 'f''(\eta)', 'f''''(\eta)'};
legend(Legend1, 'FontSize', 14);
function d=eqsol(x)
options=odeset('RelTol',1e-8,'AbsTol',[1e-8 1e-8]);
[a,b]=ode45(@eqsoll,[0,20],[0 x],options);
c=length(a);
d=b(c,1)-1;
figure(1)
plot(a,b(:,1))
xlabel('Temperature', 'FontSize', 15);
ylabel('\eta', 'FontSize', 20);
grid on
hold on;
end
end
function d_y=eqsoll(t,y)
d_y=zeros(2,1);
d_y(1)=y(2);
d_y(2)=-(0.71*18.2792*y(2)/2);
end