I am trying to code a solution to blasius eq using Runge kutta 4, help please.

clear all;

clc;

% 3 First order ODE´S from Blasius Eq

% dF/deta = G

% dG/deta = H

% dH/deta = -0.5*F*H

fF=@(eta,G) G;

fG=@(eta,H) H;

fH=@(eta,F,H) -0.5*F*H;

%initial conditions

F0 = 0;

G0 = 0;

H0 = 0; %Inital Guess for H0

% Step size and Eta max

h=0.0001;

eta=10;

N=ceil(eta/h);

%Update loop

for i=1:N

eta(i+1)=eta(i)+h;

% Runge-Kutta 4

k1F=fF(eta(i) ,F(i) ,G(i) ,H(i));

k1G=fG(eta(i) ,F(i) ,G(i) ,H(i));

k1H=fH(eta(i) ,F(i) ,G(i) ,H(i));

k2F=fF(eta(i)+h/2,F(i)+h/2*k1F,G(i)+h/2*k1G,H(i)+h/2*k1H);

k2G=fG(eta(i)+h/2,F(i)+h/2*k1F,G(i)+h/2*k1G,H(i)+h/2*k1H);

k2H=fH(eta(i)+h/2,F(i)+h/2*k1F,G(i)+h/2*k1G,H(i)+h/2*k1H);

k3F=fF(eta(i)+h/2,F(i)+h/2*k2F,G(i)+h/2*k2G,H(i)+h/2*k2H);

k3G=fG(eta(i)+h/2,F(i)+h/2*k2F,G(i)+h/2*k2G,H(i)+h/2*k2H);

k3H=fH(eta(i)+h/2,F(i)+h/2*k2F,G(i)+h/2*k2G,H(i)+h/2*k2H);

k4F=fF(eta(i)+h ,F(i)+h *k3F,G(i)+h *k3G,H(i)+h *k3H);

k4G=fG(eta(i)+h ,F(i)+h *k3F,G(i)+h *k3G,H(i)+h *k3H);

k4H=fH(eta(i)+h ,F(i)+h *k3F,G(i)+h *k3G,H(i)+h *k3H);

F(i+1)=F(i)+(h/6)*(k1F + 2*k2F + 2*k3F + k4F);

G(i+1)=G(i)+(h/6)*(k1G + 2*k2G + 2*k1G + k4G);

H(i+1)=H(i)+(h/6)*(k1G + 2*k2G + 2*k1G + k4G);

end

%Plot solution

figure(1); clf(1)

plot(eta,G)

Answer 1

[Torsten]

clear all;
clc;
% 3 First order ODE´S from Blasius Eq
% dF/deta = G
% dG/deta = H
% dH/deta = -0.5*F*H
fF=@(eta,F,G,H)          G;
fG=@(eta,F,G,H)          H;
fH=@(eta,F,G,H)   -0.5*F*H;
%initial conditions
F0 = 0;
G0 = 0;
H0 = 0; %Inital Guess for H0
F(1) = F0;
G(1) = G0;
H(1) = H0;
% Step size and Eta max
h=0.0001;
eta=10;
N=ceil(eta/h);
%Update loop
for i=1:N
    eta(i+1)=eta(i)+h;
    % Runge-Kutta 4
    k1F=fF(eta(i)    ,F(i)        ,G(i)        ,H(i));
    k1G=fG(eta(i)    ,F(i)        ,G(i)        ,H(i));
    k1H=fH(eta(i)    ,F(i)        ,G(i)        ,H(i));
    k2F=fF(eta(i)+h/2,F(i)+h/2*k1F,G(i)+h/2*k1G,H(i)+h/2*k1H);
    k2G=fG(eta(i)+h/2,F(i)+h/2*k1F,G(i)+h/2*k1G,H(i)+h/2*k1H);
    k2H=fH(eta(i)+h/2,F(i)+h/2*k1F,G(i)+h/2*k1G,H(i)+h/2*k1H);
    k3F=fF(eta(i)+h/2,F(i)+h/2*k2F,G(i)+h/2*k2G,H(i)+h/2*k2H);
    k3G=fG(eta(i)+h/2,F(i)+h/2*k2F,G(i)+h/2*k2G,H(i)+h/2*k2H);
    k3H=fH(eta(i)+h/2,F(i)+h/2*k2F,G(i)+h/2*k2G,H(i)+h/2*k2H);
    k4F=fF(eta(i)+h  ,F(i)+h  *k3F,G(i)+h  *k3G,H(i)+h  *k3H);
    k4G=fG(eta(i)+h  ,F(i)+h  *k3F,G(i)+h  *k3G,H(i)+h  *k3H);
    k4H=fH(eta(i)+h  ,F(i)+h  *k3F,G(i)+h  *k3G,H(i)+h  *k3H);
    F(i+1)=F(i)+(h/6)*(k1F + 2*k2F + 2*k3F + k4F);
    G(i+1)=G(i)+(h/6)*(k1G + 2*k2G + 2*k3G + k4G);
    H(i+1)=H(i)+(h/6)*(k1H + 2*k2H + 2*k3H + k4H);
end
%Plot solution
figure(1); clf(1)
plot(eta,G)

Answer 2

[VBBV]

@Guillermo, The anonymous functions, F ,G, H defined for the blasius flow need to applied in the same manner when RK4 method is implemented

clear all;
clc;
% 3 First order ODE´S from Blasius Eq
% dF/deta = G
% dG/deta = H
% dH/deta = -0.5*F*H

%initial conditions
F(1) = 0.01;
G(1) = 0.01;
H(1) = 0.1; %Inital Guess for H0

fF=@(eta,G)          G;
fG=@(eta,H)          H;
fH=@(eta,F,H)   -0.5*F*H;

% Step size and Eta max
h=0.0001;
eta=10;
N=ceil(eta/h);
%Update loop
for i=1:N
    eta(i+1)=eta(i)+h;
    % Runge-Kutta 4
    k1F=fF(eta(i),G(i));
    k1G=fG(eta(i),H(i));
    k1H=fH(eta(i),F(i),H(i));
    
    k2F=fF(eta(i)+h/2,G(i)+h/2*k1G);
    k2G=fG(eta(i)+h/2,H(i)+h/2*k1H);
    k2H=fH(eta(i)+h/2,F(i)+h/2*k1F,H(i)+h/2*k1H);
    
    k3F=fF(eta(i)+h/2,G(i)+h/2*k2G);
    k3G=fG(eta(i)+h/2,H(i)+h/2*k2H);    
    k3H=fH(eta(i)+h/2,F(i)+h/2*k2F,H(i)+h/2*k2H);
    
    k4F=fF(eta(i)+h,G(i)+h*k3G);
    k4G=fG(eta(i)+h,H(i)+h*k3H);
    k4H=fH(eta(i)+h,F(i)+h*k3F,H(i)+h*k3H);
    
    F(i+1)=F(i)+(h/6)*(k1F + 2*k2F + 2*k3F + k4F);
    G(i+1)=G(i)+(h/6)*(k1G + 2*k2G + 2*k1G + k4G);
    H(i+1)=H(i)+(h/6)*(k1H + 2*k2H + 2*k1H + k4H);
end
%Plot solution
hold on
subplot(311);plot(eta,F); subplot(312); plot(eta,G); subplot(313);plot(eta,H);