while solving a couple ordinary differential equation, I'm getting an error of Unable to perform assignment because the left and right sides have a different number of element

output:-

Unable to perform assignment because the left and right sides have a different number of element.

Error in laser2 (line 52)

A1(i+1) = A1(i) + (h/6)*(k1+(2*(k2+k3))+k4);

input:-

format long
a = 0;
b = 30;
h = 0.1;
o = 3.5;  
tc = 3.0;  
tf = 2.3;  
a1 = 0.1;
a2 = 0.1;
P1 = 0.2;  
P2 = 0.2;
k= 0.9;
V = 1;
A1(1) = 10E-4; 
A2(1) = 10E-5;
G1(1) = 10E-6;
G2(1) = 10E-6;
O(1) = 0;
n = (b-a)/h ;
t = ((a) + (0:n)*h).*10E-6;
func1 = @(G1,A1,A2,O) (1/tc).*(G1- a1).*(A1./V) +(k./tc).*(A2./V).*cos(O);           
func3 = @(G2,A2,A1,O)  (1/tc).*(G2- a2).*(A2./V) +(k./tc).*(A1./V).*cos(O);            
func2 = @(G1)    (1/tf).*(P1 - G1.*(A1.^2+1));
func4 = @(G2)   (1/tf).*(P2 - G2.*(A2.^2+1));
func5 = @(A1,A2,O)   o - (k./tc).*((A1./A2) + (A2./A1)).*sin(O);
for i = 1:n
    k1 = feval(func1,G1(i),A1(i),A2(i),O(i));
    l1 = feval(func2, G1(i));
    m1 = feval(func3, G2(i),A1(i),A2(i),O(i));
    n1 = feval(func4, G2(i));
    q1 = feval(func5,A1(i),A2(i),O(i));
    
    k2 = feval(func1, A2(i)+(m1/2)*h, A1(i)+(k1/2)*h,  G1(i)+(l1/2)*h,  O(i)+(q1/2)*h);
    l2 = feval(func2, G1(i)+(l1*h/2));
    m2 = feval(func3, A2(i)+(m1/2)*h, A1(i)+(k1/2)*h,  G2(i)+(n1/2)*h,  O(i)+(q1/2)*h);
    n2 = feval(func4, G2(i)+(n1/2)*h);
    q2 = feval(func5, A2(i)+(m1/2)*h, A1(i)+(k1/2)*h, O(i)+(q1/2)*h);
    
    k3 = feval(func1, A2(i)+(m2/2)*h, A1(i)+(k2/2)*h,  G1(i)+(l2/2)*h,  O(i)+(q2/2)*h);
    l3 = feval(func2, G1(i)+(l2*h/2));
    m3 = feval(func3, A2(i)+(m2/2)*h, A1(i)+(k2/2)*h,  G2(i)+(n2/2)*h,  O(i)+(q2/2)*h);
    n3 = feval(func4, G2(i)+(n2/2)*h);
    q3 = feval(func5, A2(i)+(m2/2)*h, A1(i)+(k2/2)*h, O(i)+(q2/2)*h );
    
    k4 = feval(func1, A2(i)+(m3*h), A1(i)+(k3*h),  G1(i)+(l3*h),  O(i)+(q3*h));
    l4 = feval(func2, G1(i)+(l3*h));
    m4 = feval(func3, A2(i)+(m3*h), A1(i)+(k3*h),  G2(i)+(n3*h),  O(i)+(q3*h));
    n4 = feval(func4, G2(i)+(n3*h));
    q4 = feval(func5, A2(i)+(m3*h), A1(i)+(k3*h), O(i)+(q3*h));
    
    A1(i+1) = A1(i) + (h/6)*(k1+(2*(k2+k3))+k4);   
    G1(i+1) = G1(i) + (h/6)*(l1+(2*(l2+l3))+l4);
    A2(i+1) = A2(i) + (h/6)*(m1+(2*(m2+m3))+m4); 
    G2(i+1) = G2(i) + (h/6)*(n1+(2*(n2+n3))+n4); 
    O(i+1)  = O(i) + (h/6)*(q1+(2*(q2+q3))+q4);
end
subplot 211
plot(t,A1)
grid on
subplot 212
plot(t,O)

it maybe due to the func2 = @(G1) (1/tf).*(P1 - G1.*(A1.^2+1));

func4 = @(G2) (1/tf).*(P2 - G2.*(A2.^2+1));

this two code, but I couldn't understand what's the problem in it.

Answer 1

[Johan]

func2 and func4 uses A1 and A2 respectively which means they are using the whole arrays, either call A1 and A2 in the func definition (see example below) or set it to A1(1) and A2(1) if you only want to use the first term.

On a side note it's considered bad practice to add 1 2 3 .... to your variable names. It makes your code hard to understand and debug there is a post that explain why here : https://fr.mathworks.com/support/search.html/answers/304528-tutorial-why-variables-should-not-be-named-dynamically-eval.html

%change some variable names so that they make sense
time_start = 0;
time_end = 30;
time_step = 0.1;
o = 3.5;  
tc = 3.0;  
tf = 2.3;  
a1 = 0.1;
a2 = 0.1;
P1 = 0.2;  
P2 = 0.2;
k= 0.9;
V = 1;
n_step = (time_end-time_start)/time_step ;
time = ((time_start) + (0:n_step)*time_step).*10E-6;
func1 = @(G1,A1,A2,O) (1/tc).*(G1- a1).*(A1./V) +(k./tc).*(A2./V).*cos(O);           
func3 = @(G2,A2,A1,O)  (1/tc).*(G2- a2).*(A2./V) +(k./tc).*(A1./V).*cos(O);            
func2 = @(G1,A1)    (1/tf).*(P1 - G1.*(A1.^2+1));
func4 = @(G2,A2)   (1/tf).*(P2 - G2.*(A2.^2+1));
func5 = @(A1,A2,O)   o - (k./tc).*((A1./A2) + (A2./A1)).*sin(O);

%Initialize arrays
A1 = zeros(n_step,1);
A2 = A1;
G1 = A1;
G2 = A1;
O = A1;

A1(1) = 10E-4; 
A2(1) = 10E-5;
G1(1) = 10E-6;
G2(1) = 10E-6;
O(1) = 0;

for i_loop = 1:n_step
    k1 = func1(G1(i_loop),A1(i_loop),A2(i_loop),O(i_loop));
    l1 = func2( G1(i_loop),A1(i_loop));
    m1 = func3( G2(i_loop),A1(i_loop),A2(i_loop),O(i_loop));
    n1 = func4( G2(i_loop),A2(i_loop));
    q1 = func5(A1(i_loop),A2(i_loop),O(i_loop));
    
    k2 = func1( A2(i_loop)+(m1/2)*time_step, A1(i_loop)+(k1/2)*time_step,  G1(i_loop)+(l1/2)*time_step,  O(i_loop)+(q1/2)*time_step);
    l2 = func2( G1(i_loop)+(l1*time_step/2), A1(i_loop)+(k1/2)*time_step);
    m2 = func3( A2(i_loop)+(m1/2)*time_step, A1(i_loop)+(k1/2)*time_step,  G2(i_loop)+(n1/2)*time_step,  O(i_loop)+(q1/2)*time_step);
    n2 = func4( G2(i_loop)+(n1/2)*time_step, A2(i_loop)+(m1/2)*time_step);
    q2 = func5( A2(i_loop)+(m1/2)*time_step, A1(i_loop)+(k1/2)*time_step, O(i_loop)+(q1/2)*time_step);
    
    k3 = func1( A2(i_loop)+(m2/2)*time_step, A1(i_loop)+(k2/2)*time_step,  G1(i_loop)+(l2/2)*time_step,  O(i_loop)+(q2/2)*time_step);
    l3 = func2( G1(i_loop)+(l2*time_step/2), A1(i_loop)+(k2/2)*time_step);
    m3 = func3( A2(i_loop)+(m2/2)*time_step, A1(i_loop)+(k2/2)*time_step,  G2(i_loop)+(n2/2)*time_step,  O(i_loop)+(q2/2)*time_step);
    n3 = func4( G2(i_loop)+(n2/2)*time_step, A2(i_loop)+(m2/2)*time_step);
    q3 = func5( A2(i_loop)+(m2/2)*time_step, A1(i_loop)+(k2/2)*time_step, O(i_loop)+(q2/2)*time_step );
    
    k4 = func1( A2(i_loop)+(m3*time_step), A1(i_loop)+(k3*time_step),  G1(i_loop)+(l3*time_step),  O(i_loop)+(q3*time_step));
    l4 = func2( G1(i_loop)+(l3*time_step), A1(i_loop)+(k3*time_step));
    m4 = func3( A2(i_loop)+(m3*time_step), A1(i_loop)+(k3*time_step),  G2(i_loop)+(n3*time_step),  O(i_loop)+(q3*time_step));
    n4 = func4( G2(i_loop)+(n3*time_step),A2(i_loop)+(m3*time_step));
    q4 = func5( A2(i_loop)+(m3*time_step), A1(i_loop)+(k3*time_step), O(i_loop)+(q3*time_step));
    
    A1(i_loop+1) = A1(i_loop) + (time_step/6)*(k1+(2*(k2+k3))+k4);   
    G1(i_loop+1) = G1(i_loop) + (time_step/6)*(l1+(2*(l2+l3))+l4);
    A2(i_loop+1) = A2(i_loop) + (time_step/6)*(m1+(2*(m2+m3))+m4); 
    G2(i_loop+1) = G2(i_loop) + (time_step/6)*(n1+(2*(n2+n3))+n4); 
    O(i_loop+1)  = O(i_loop) + (time_step/6)*(q1+(2*(q2+q3))+q4);
end
subplot 211
plot(time,A1)
grid on
subplot 212
plot(time,O)

Answer 2

[KSSV]

This function:

func2 = @(G1)    (1/tf).*(P1 - G1.*(A1.^2+1));

should be:

func2 = @(G1,A1)    (1/tf).*(P1 - G1.*(A1.^2+1));

Answer 3

[Walter Roberson]

The posted code executes for me

format long
a = 0;
b = 30;
h = 0.1;
o = 3.5;  
tc = 3.0;  
tf = 2.3;  
a1 = 0.1;
a2 = 0.1;
P1 = 0.2;  
P2 = 0.2;
k= 0.9;
V = 1;
A1(1) = 10E-4; 
A2(1) = 10E-5;
G1(1) = 10E-6;
G2(1) = 10E-6;
O(1) = 0;
n = (b-a)/h ;
t = ((a) + (0:n)*h).*10E-6;
func1 = @(G1,A1,A2,O) (1/tc).*(G1- a1).*(A1./V) +(k./tc).*(A2./V).*cos(O);           
func3 = @(G2,A2,A1,O)  (1/tc).*(G2- a2).*(A2./V) +(k./tc).*(A1./V).*cos(O);            
func2 = @(G1)    (1/tf).*(P1 - G1.*(A1.^2+1));
func4 = @(G2)   (1/tf).*(P2 - G2.*(A2.^2+1));
func5 = @(A1,A2,O)   o - (k./tc).*((A1./A2) + (A2./A1)).*sin(O);
for i = 1:n
    k1 = feval(func1,G1(i),A1(i),A2(i),O(i));
    l1 = feval(func2, G1(i));
    m1 = feval(func3, G2(i),A1(i),A2(i),O(i));
    n1 = feval(func4, G2(i));
    q1 = feval(func5,A1(i),A2(i),O(i));
    
    k2 = feval(func1, A2(i)+(m1/2)*h, A1(i)+(k1/2)*h,  G1(i)+(l1/2)*h,  O(i)+(q1/2)*h);
    l2 = feval(func2, G1(i)+(l1*h/2));
    m2 = feval(func3, A2(i)+(m1/2)*h, A1(i)+(k1/2)*h,  G2(i)+(n1/2)*h,  O(i)+(q1/2)*h);
    n2 = feval(func4, G2(i)+(n1/2)*h);
    q2 = feval(func5, A2(i)+(m1/2)*h, A1(i)+(k1/2)*h, O(i)+(q1/2)*h);
    
    k3 = feval(func1, A2(i)+(m2/2)*h, A1(i)+(k2/2)*h,  G1(i)+(l2/2)*h,  O(i)+(q2/2)*h);
    l3 = feval(func2, G1(i)+(l2*h/2));
    m3 = feval(func3, A2(i)+(m2/2)*h, A1(i)+(k2/2)*h,  G2(i)+(n2/2)*h,  O(i)+(q2/2)*h);
    n3 = feval(func4, G2(i)+(n2/2)*h);
    q3 = feval(func5, A2(i)+(m2/2)*h, A1(i)+(k2/2)*h, O(i)+(q2/2)*h );
    
    k4 = feval(func1, A2(i)+(m3*h), A1(i)+(k3*h),  G1(i)+(l3*h),  O(i)+(q3*h));
    l4 = feval(func2, G1(i)+(l3*h));
    m4 = feval(func3, A2(i)+(m3*h), A1(i)+(k3*h),  G2(i)+(n3*h),  O(i)+(q3*h));
    n4 = feval(func4, G2(i)+(n3*h));
    q4 = feval(func5, A2(i)+(m3*h), A1(i)+(k3*h), O(i)+(q3*h));
    
    A1(i+1) = A1(i) + (h/6)*(k1+(2*(k2+k3))+k4);   
    G1(i+1) = G1(i) + (h/6)*(l1+(2*(l2+l3))+l4);
    A2(i+1) = A2(i) + (h/6)*(m1+(2*(m2+m3))+m4); 
    G2(i+1) = G2(i) + (h/6)*(n1+(2*(n2+n3))+n4); 
    O(i+1)  = O(i) + (h/6)*(q1+(2*(q2+q3))+q4);
end
subplot 211
plot(t,A1)
grid on
subplot 212
plot(t,O)