% system of first order differential equation
% y1' = y2;
% y2' = y3;
% y3' = (1-phi).^2.5*[(1-phi)+phi*(rows/rowp)*(y1*y3-y1.^2)-M*y2+lamda*y4*((1-phi)+phi*(rowbetas/rowbetaf))]
% y4' = y5;
% y5' = -(3*N*pr*kf*[(1-phi)+phi*(rowcps/rowcpf)]/(3*N+4)*knf)*y1*y5;
%with boundary condition y1(0) = 0, y2(0)=1, y3(0)=0.078, y4(0) = 1, y5(0)=?
phi = 0.01;
lamda = 1;
M = 2;
Ec = 2;
N = 1.5;
pr = 6.2;
%copper
rows =8933 ;
ks = 401;
cps =385 ;
rowcps=8933*385;
rowbetas = 8933*1.67*10.^-5;
betas = 1.67*10.^-5;
%alumina
% rows =3970 ;
% ks =40 ;
% cps =765;
%betas = 0.85*10.^-5;
%water
rowf = 997.1;
kf = 0.613;
cpf = 4179;
betaf = 21*10.^-5;
rowbetaf = 997.1*21*10.^-5;
rowcpf = 997.1*4179;
% constants
A1 = (1-phi)+phi*(rows/rowf);
A2 = (1-phi)+phi*(rowbetas/rowbetaf);
A3 = (1-phi)+phi*(rowcps/rowcpf);
A = (1-phi).^2.5;
A4 = -(3*N*pr*kf/(3*N+4)*knf);
knf = kf*(ks+2*kf-2*phi*(kf-ks)/ks+2*kf+phi*(kf-ks));
%% solve ODE-IVP using RK-4 Method
tic;
% range of independent variable η
eta_0 = 0;
eta_max = 5;
% Initial conditions
f0 = 0;
g0 = 1;
h0 = -0.74878322
l0 = 1;
m0 = 0.2356
% Step value Δη
d_eta = 0.001;
N = (eta_max - eta_0)/ d_eta;
err= 1;
%% Initializing solution
eta = eta_0 : d_eta : eta_max; % X coordinate for plot function
while err >10^-9
F = zeros(N+1,3); % f = F(N+1,1) ; g = F(N+1,2) ; h = F(N+1,3) ;
F(1,:) = [f0 g0 h0];
G = zeros(N+1,2);
G(1,:) = [l0 m0]; %l = G(N+1,1); m = G(N+1,2); j let us choose artbitrary for two values
%% Trail using RK-4 Method
for i = 1:N
k11(i,:) = d_eta*[G(i,2) A4*A3*F(i,1)*G(i,2)];
k1(i,:) = d_eta* [F(i,2) F(i,3) A*(A1*(F(i,1)*F(i,3)-F(i,2).^2)-M*F(i,2)+lamda*A2*G(i,1))] ;
k12(i,:) = d_eta*[G(i,2)+k11(i,1)/2 A4*A3*F(i,1)+k1(i,1)/2*G(i,2)+k11(i,2)/2];
k2(i,:) = d_eta* [F(i,2)+k1(i,2)/2 F(i,3)+k1(i,3)/2 A*(A1*(F(i,1)+k1(i,1)/2*F(i,3)+k1(i,3)/2-F(i,2).^2+k1(i,2)/2)-M*F(i,2)+k1(i,2)/2+lamda*A2*G(i,1)+k11(i,1)) ];
k13(i,:) = d_eta*[G(i,2)+k12(i,1)/2 A4*A3*F(i,1)+k2(i,1)/2*G(i,2)+k12(i,2)/2];
k3(i,:) = d_eta* [F(i,2)+k2(i,2)/2 F(i,3)+k2(i,3)/2 A*(A1*(F(i,1)+k2(i,1)/2*F(i,3)+k2(i,3)/2-F(i,2).^2+k2(i,2)/2)-M*F(i,2)+k2(i,2)/2+lamda*A2*G(i,1)+k12(i,1)/2) ];
k14(i,:) = d_eta*[G(i,2)+k13(i,1) A4*A3*F(i,1)+k3(i,1)*G(i,2)+k13(i,2)];
k4(i,:) = d_eta* [F(i,2)+k3(i,2) F(i,3)+k3(i,3) A*(A1*(F(i,1)+k3(i,1)*F(i,3)+k3(i,3)-F(i,2).^2+k3(i,2))-M*F(i,2)+k3(i,2)+lamda*A2*G(i,1)+k13(i,1))];
G(i+1,:) = G(i,:)+(1/6)*(k11(i,:)+2*k12(i,:)+2*k13(i,:)+k14(i,:));
F(i+1,:) = F(i,:)+(1/6)*(k1(i,:)+2*k2(i,:)+2*k3(i,:)+k4(i,:));
end
end
%% plotting
f1=figure;
f1.Units = 'normalized';
f1.Position = [0.1 0.1 0.8 0.6];
%plot(eta',F,'LineWidth',1);
plot(eta',G,'Linewidth',1);
