Unable to perform assignment because the size of the left side is 1-by-3 and the size of the right side is 1-by-2.

% 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);

Answer 1

[Walter Roberson]

d_eta = 0.001;
N = (eta_max - eta_0)/ d_eta;

You use N as a size, so you are assuming that it is an integer. However, 0.001 is not exactly representatable in IEEE 754 double precision floating point, so dividing by 0.001 is not certain to give you back an integer, and is not exactly the same as multiplying by 1000.