I’m using Runge-Kutta method to solve a coupled system of EDO.
The function y1 is valid until a certain time (here is the "ts2" parameter). After this time, y8 is valid.
I want to plot y1 from 0 to ts2 (which happens without problem) and y8 from ts2 until the end.
But, I get the error "Arrays have incompatible sizes for this operation." during the second part.
%Parameters
or =0.2;
oc = 0.2;

gamma = 0.01;
gam=0.01;
Gamma = 1;
Gamma1 = Gamma;
Gamma2 = Gamma;
Gamma3 = 9/10;
Gamma5 = 1/10;

sig13ts2 = 5.2;
sig57ts = 2.6;
ts2=10;

f1 = @(t,y1,y2) (-Gamma/2-gam/2)*y1 -i*or*y2;
f2 = @(t,y1,y2) -i*or*y1+(-gam/2)*y2;

h=0.0001;
t(1) = 0;
y1(1) = 0;
y2(1) = sig57ts;
for j=1:1100000
t(j+1)=t(j)+h;
k1y1 = h*f1(t(j), y1(j), y2(j));
k1y2 = h*f2(t(j), y1(j), y2(j));

k2y1 = h*f1(t(j)+h/2,y1(j)+k1y1/2,y2(j)+k1y2/2);
k2y2 = h*f2(t(j)+h/2,y1(j)+k1y1/2,y2(j)+k1y2/2);

k3y1 = h*f1(t(j)+h/2,y1(j)+k2y1/2,y2(j)+k2y2/2);
k3y2 = h*f2(t(j)+h/2,y1(j)+k2y1/2,y2(j)+k2y2/2);

k4y1 = h*f1(t(j)+h, y1(j)+k3y1, y2(j)+k3y2);
k4y2 = h*f2(t(j)+h, y1(j)+k3y1, y2(j)+k3y2);

y1(j+1)=y1(j) + (k1y1 + 2*k2y1 + 2*k3y1 + k4y1)/6;
y2(j+1)=y2(j) + (k1y2 + 2*k2y2 + 2*k3y2 + k4y2)/6;
end
plot(t(1:ts2/h),abs(y1(1:ts2/h)).^2)
xlim([0 60])
hold on

%% Ligando C apos ts2
f3 = @(t,y3,y4,y5,y6,y7,y8) -(+Gamma2/2+gam/2)*y3+1i*oc*y4-1i*oc*y5;
f4 = @(t,y3,y4,y5,y6,y7,y8) 1i*oc*y3-(+gamma/2)*y4+Gamma3*y5-1i*oc*y6;
f5 = @(t,y3,y4,y5,y6,y7,y8) -1i*oc*y3-(Gamma1)*y5+1i*oc*y6;
f6 = @(t,y3,y4,y5,y6,y7,y8) -1i*oc*y4+1i*oc*y5+(-Gamma1/2-gamma/2)*y6;
f7 = @(t,y3,y4,y5,y6,y7,y8) -1i*or*y8+(-Gamma2/2-gam/2)*y7;
f8 = @(t,y3,y4,y5,y6,y7,y8) -1i*or*y7+Gamma5*y5-(gam/2)*y8;

t(1) = ts2;
y3(1) = 0;
y4(1) = sig13ts2;
y5(1) = 0;
y6(1) = 0;
y7(1) = y1(ts2/h);
y8(1) = y2(ts2/h);

for j=1:1000000
t(j+1)=t(j)+h;
k1y3 = h*f3(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y4 = h*f4(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y5 = h*f5(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y6 = h*f6(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y7 = h*f7(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y8 = h*f8(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));

k2y3 = h*f3(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y4 = h*f4(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y5 = h*f5(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y6 = h*f6(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y7 = h*f7(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y8 = h*f8(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);

k3y3 = h*f3(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y4 = h*f4(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y5 = h*f5(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y6 = h*f6(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y7 = h*f7(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y8 = h*f8(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);

k4y3 = h*f3(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y4 = h*f4(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y5 = h*f5(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y6 = h*f6(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y7 = h*f7(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y8 = h*f8(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);

y3(j+1)=y3(j) + (k1y3 + 2*k2y3 + 2*k3y3 + k4y3)/6;
y4(j+1)=y4(j) + (k1y4 + 2*k2y4 + 2*k3y4 + k4y4)/6;
y5(j+1)=y5(j) + (k1y5 + 2*k2y5 + 2*k3y5 + k4y5)/6;
y6(j+1)=y6(j) + (k1y6 + 2*k2y6 + 2*k3y6 + k4y6)/6;
y7(j+1)=y7(j) + (k1y7 + 2*k2y7 + 2*k3y7 + k4y7)/6;
y8(j+1)=y8(j) + (k1y8 + 2*k2y8 + 2*k3y8 + k4y8)/6;
end

plot(t,abs(y8).^2+abs(y1(ts2/h:end)).^2)I’m using Runge-Kutta method to solve a coupled system of EDO.
The function y1 is valid until a certain time (here is the "ts2" parameter). After this time, y8 is valid.
I want to plot y1 from 0 to ts2 (which happens without problem) and y8 from ts2 until the end.
But, I get the error "Arrays have incompatible sizes for this operation." during the second part.
%Parameters
or =0.2;
oc = 0.2;

gamma = 0.01;
gam=0.01;
Gamma = 1;
Gamma1 = Gamma;
Gamma2 = Gamma;
Gamma3 = 9/10;
Gamma5 = 1/10;

sig13ts2 = 5.2;
sig57ts = 2.6;
ts2=10;

f1 = @(t,y1,y2) (-Gamma/2-gam/2)*y1 -i*or*y2;
f2 = @(t,y1,y2) -i*or*y1+(-gam/2)*y2;

h=0.0001;
t(1) = 0;
y1(1) = 0;
y2(1) = sig57ts;
for j=1:1100000
t(j+1)=t(j)+h;
k1y1 = h*f1(t(j), y1(j), y2(j));
k1y2 = h*f2(t(j), y1(j), y2(j));

k2y1 = h*f1(t(j)+h/2,y1(j)+k1y1/2,y2(j)+k1y2/2);
k2y2 = h*f2(t(j)+h/2,y1(j)+k1y1/2,y2(j)+k1y2/2);

k3y1 = h*f1(t(j)+h/2,y1(j)+k2y1/2,y2(j)+k2y2/2);
k3y2 = h*f2(t(j)+h/2,y1(j)+k2y1/2,y2(j)+k2y2/2);

k4y1 = h*f1(t(j)+h, y1(j)+k3y1, y2(j)+k3y2);
k4y2 = h*f2(t(j)+h, y1(j)+k3y1, y2(j)+k3y2);

y1(j+1)=y1(j) + (k1y1 + 2*k2y1 + 2*k3y1 + k4y1)/6;
y2(j+1)=y2(j) + (k1y2 + 2*k2y2 + 2*k3y2 + k4y2)/6;
end
plot(t(1:ts2/h),abs(y1(1:ts2/h)).^2)
xlim([0 60])
hold on

%% Ligando C apos ts2
f3 = @(t,y3,y4,y5,y6,y7,y8) -(+Gamma2/2+gam/2)*y3+1i*oc*y4-1i*oc*y5;
f4 = @(t,y3,y4,y5,y6,y7,y8) 1i*oc*y3-(+gamma/2)*y4+Gamma3*y5-1i*oc*y6;
f5 = @(t,y3,y4,y5,y6,y7,y8) -1i*oc*y3-(Gamma1)*y5+1i*oc*y6;
f6 = @(t,y3,y4,y5,y6,y7,y8) -1i*oc*y4+1i*oc*y5+(-Gamma1/2-gamma/2)*y6;
f7 = @(t,y3,y4,y5,y6,y7,y8) -1i*or*y8+(-Gamma2/2-gam/2)*y7;
f8 = @(t,y3,y4,y5,y6,y7,y8) -1i*or*y7+Gamma5*y5-(gam/2)*y8;

t(1) = ts2;
y3(1) = 0;
y4(1) = sig13ts2;
y5(1) = 0;
y6(1) = 0;
y7(1) = y1(ts2/h);
y8(1) = y2(ts2/h);

for j=1:1000000
t(j+1)=t(j)+h;
k1y3 = h*f3(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y4 = h*f4(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y5 = h*f5(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y6 = h*f6(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y7 = h*f7(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y8 = h*f8(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));

k2y3 = h*f3(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y4 = h*f4(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y5 = h*f5(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y6 = h*f6(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y7 = h*f7(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y8 = h*f8(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);

k3y3 = h*f3(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y4 = h*f4(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y5 = h*f5(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y6 = h*f6(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y7 = h*f7(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y8 = h*f8(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);

k4y3 = h*f3(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y4 = h*f4(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y5 = h*f5(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y6 = h*f6(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y7 = h*f7(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y8 = h*f8(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);

y3(j+1)=y3(j) + (k1y3 + 2*k2y3 + 2*k3y3 + k4y3)/6;
y4(j+1)=y4(j) + (k1y4 + 2*k2y4 + 2*k3y4 + k4y4)/6;
y5(j+1)=y5(j) + (k1y5 + 2*k2y5 + 2*k3y5 + k4y5)/6;
y6(j+1)=y6(j) + (k1y6 + 2*k2y6 + 2*k3y6 + k4y6)/6;
y7(j+1)=y7(j) + (k1y7 + 2*k2y7 + 2*k3y7 + k4y7)/6;
y8(j+1)=y8(j) + (k1y8 + 2*k2y8 + 2*k3y8 + k4y8)/6;
end

plot(t,abs(y8).^2+abs(y1(ts2/h:end)).^2) I’m using Runge-Kutta method to solve a coupled system of EDO.
The function y1 is valid until a certain time (here is the "ts2" parameter). After this time, y8 is valid.
I want to plot y1 from 0 to ts2 (which happens without problem) and y8 from ts2 until the end.
But, I get the error "Arrays have incompatible sizes for this operation." during the second part.
%Parameters
or =0.2;
oc = 0.2;

gamma = 0.01;
gam=0.01;
Gamma = 1;
Gamma1 = Gamma;
Gamma2 = Gamma;
Gamma3 = 9/10;
Gamma5 = 1/10;

sig13ts2 = 5.2;
sig57ts = 2.6;
ts2=10;

f1 = @(t,y1,y2) (-Gamma/2-gam/2)*y1 -i*or*y2;
f2 = @(t,y1,y2) -i*or*y1+(-gam/2)*y2;

h=0.0001;
t(1) = 0;
y1(1) = 0;
y2(1) = sig57ts;
for j=1:1100000
t(j+1)=t(j)+h;
k1y1 = h*f1(t(j), y1(j), y2(j));
k1y2 = h*f2(t(j), y1(j), y2(j));

k2y1 = h*f1(t(j)+h/2,y1(j)+k1y1/2,y2(j)+k1y2/2);
k2y2 = h*f2(t(j)+h/2,y1(j)+k1y1/2,y2(j)+k1y2/2);

k3y1 = h*f1(t(j)+h/2,y1(j)+k2y1/2,y2(j)+k2y2/2);
k3y2 = h*f2(t(j)+h/2,y1(j)+k2y1/2,y2(j)+k2y2/2);

k4y1 = h*f1(t(j)+h, y1(j)+k3y1, y2(j)+k3y2);
k4y2 = h*f2(t(j)+h, y1(j)+k3y1, y2(j)+k3y2);

y1(j+1)=y1(j) + (k1y1 + 2*k2y1 + 2*k3y1 + k4y1)/6;
y2(j+1)=y2(j) + (k1y2 + 2*k2y2 + 2*k3y2 + k4y2)/6;
end
plot(t(1:ts2/h),abs(y1(1:ts2/h)).^2)
xlim([0 60])
hold on

%% Ligando C apos ts2
f3 = @(t,y3,y4,y5,y6,y7,y8) -(+Gamma2/2+gam/2)*y3+1i*oc*y4-1i*oc*y5;
f4 = @(t,y3,y4,y5,y6,y7,y8) 1i*oc*y3-(+gamma/2)*y4+Gamma3*y5-1i*oc*y6;
f5 = @(t,y3,y4,y5,y6,y7,y8) -1i*oc*y3-(Gamma1)*y5+1i*oc*y6;
f6 = @(t,y3,y4,y5,y6,y7,y8) -1i*oc*y4+1i*oc*y5+(-Gamma1/2-gamma/2)*y6;
f7 = @(t,y3,y4,y5,y6,y7,y8) -1i*or*y8+(-Gamma2/2-gam/2)*y7;
f8 = @(t,y3,y4,y5,y6,y7,y8) -1i*or*y7+Gamma5*y5-(gam/2)*y8;

t(1) = ts2;
y3(1) = 0;
y4(1) = sig13ts2;
y5(1) = 0;
y6(1) = 0;
y7(1) = y1(ts2/h);
y8(1) = y2(ts2/h);

for j=1:1000000
t(j+1)=t(j)+h;
k1y3 = h*f3(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y4 = h*f4(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y5 = h*f5(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y6 = h*f6(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y7 = h*f7(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));
k1y8 = h*f8(t(j), y3(j), y4(j), y5(j), y6(j), y7(j), y8(j));

k2y3 = h*f3(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y4 = h*f4(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y5 = h*f5(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y6 = h*f6(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y7 = h*f7(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);
k2y8 = h*f8(t(j)+h/2,y3(j)+k1y3/2,y4(j)+k1y4/2,y5(j)+k1y5/2,y6(j)+k1y6/2,y7(j)++k1y7/2,y8(j)+k1y8/2);

k3y3 = h*f3(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y4 = h*f4(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y5 = h*f5(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y6 = h*f6(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y7 = h*f7(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);
k3y8 = h*f8(t(j)+h/2,y3(j)+k2y3/2,y4(j)+k2y4/2,y5(j)+k2y5/2,y6(j)+k2y6/2,y7(j)+k2y7/2,y8(j)+k2y8/2);

k4y3 = h*f3(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y4 = h*f4(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y5 = h*f5(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y6 = h*f6(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y7 = h*f7(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);
k4y8 = h*f8(t(j)+h, y3(j)+k3y3, y4(j)+k3y4,y5(j)+k3y5,y6(j)+k3y6,y7(j)+k3y7,y8(j)+k3y8);

y3(j+1)=y3(j) + (k1y3 + 2*k2y3 + 2*k3y3 + k4y3)/6;
y4(j+1)=y4(j) + (k1y4 + 2*k2y4 + 2*k3y4 + k4y4)/6;
y5(j+1)=y5(j) + (k1y5 + 2*k2y5 + 2*k3y5 + k4y5)/6;
y6(j+1)=y6(j) + (k1y6 + 2*k2y6 + 2*k3y6 + k4y6)/6;
y7(j+1)=y7(j) + (k1y7 + 2*k2y7 + 2*k3y7 + k4y7)/6;
y8(j+1)=y8(j) + (k1y8 + 2*k2y8 + 2*k3y8 + k4y8)/6;
end

plot(t,abs(y8).^2+abs(y1(ts2/h:end)).^2) vector, plot, numerical integration MATLAB Answers — New Questions

​