Plotting sum of two vectors

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)

Arrays have incompatible sizes for this operation.

Answer 1

[Torsten]

%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+1),abs(y1(1:ts2/h+1)).^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 = [];
t(1) = ts2;
y3(1) = 0;
y4(1) = sig13ts2;
y5(1) = 0;
y6(1) = 0;
y7(1) = y1(ts2/h+1);
y8(1) = y2(ts2/h+1);

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+1:end)).^2)
hold off

●

[Eduardo]

Thank you for answering, it was really helpful.

Answer 2

[Stephen23]

Use numbered variable names if you really enjoy writing lots of code.

Otherwise do something more like this:

%% 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;
h = 0.0001;

%% Phase 1: 2-variable system (t = 0 to ts2)
% dy/dt = A*y where y = [y1; y2]
A1 = [-Gamma/2-gam/2, -1i*or; -1i*or, -gam/2];

% Initial conditions
t1 = 0:h:ts2;
n1 = length(t1);
y = zeros(2,n1);
y(:,1) = [0; sig57ts];

% RK4 integration for phase 1
for j = 1:n1-1
    k1 = h * A1 * y(:,j);
    k2 = h * A1 * (y(:,j) + k1/2);
    k3 = h * A1 * (y(:,j) + k2/2);
    k4 = h * A1 * (y(:,j) + k3);
    y(:,j+1) = y(:,j) + (k1 + 2*k2 + 2*k3 + k4)/6;
end

% Plot phase 1
plot(t1, abs(y(1,:)).^2)
xlim([0,60])
hold on

%% Phase 2: 6-variable system (t = ts2 onwards)
% System matrix for [y3; y4; y5; y6; y7; y8]
A2 = [-Gamma2/2-gam/2,   1i*oc,   -1i*oc,                 0,               0,       0;...
                1i*oc, -gamma/2,  Gamma3,            -1i*oc,               0,       0;...
               -1i*oc,        0, -Gamma1,             1i*oc,               0,       0;...
                    0,   -1i*oc,   1i*oc, -Gamma1/2-gamma/2,               0,       0;...
                    0,        0,       0,                 0, -Gamma2/2-gam/2,  -1i*or;...
                    0,        0,  Gamma5,                 0,           -1i*or, -gam/2];

% Initial conditions for phase 2
n2 = 1000000;
t2 = zeros(1,n2+1);
t2(1) = ts2;
z = zeros(6, n2+1);
z(:,1) = [0; sig13ts2; 0; 0; y(1,end); y(2,end)];

% RK4 integration for phase 2
for j = 1:n2
    t2(j+1) = t2(j) + h;
    k1 = h * A2 * z(:,j);
    k2 = h * A2 * (z(:,j) + k1/2);
    k3 = h * A2 * (z(:,j) + k2/2);
    k4 = h * A2 * (z(:,j) + k3);
    z(:,j+1) = z(:,j) + (k1 + 2*k2 + 2*k3 + k4)/6;
end

plot(t2, abs(z(5,:)).^2 + abs(z(6,:)).^2)