Vector length is not working for my ODE graph code, anyone know how to resolve this issue? Many thanks

clc; clear all; close all;

opts = odeset('RelTol',1e-6,'AbsTol',[1e-9 1e-9]);    % setting ODE solver tolerances

%% setting parameters and initial conditions

r= (0.1:0.1:1.5) ;
zeta= 0.35 ;
epsilon = -1/6 ;
Y0= [0.50, 0] ;
tspan= [0 1000] ;
f = 0.45 ;

opts = odeset('RelTol',1e-6,'AbsTol',[1e-9 1e-9]);    % setting ODE solver tolerances

%% ODE solver
for i = (1:1:15) ;
    
    [t,y] = ode45(@(t,y) odefcn_script(t,y,f,r(i),zeta), tspan, Y0);
    
    
    harmonic_equation = @(A) f^2-(-r(i)^2*A + A +((epsilon*3*A^3)/4)^2)-((2*A*zeta*r(i))^2);
    
    amplitude(:,i) = fzero(harmonic_equation, 1);   %solving harmonic balace equation
    
    theta1 = amplitude(i)*sin(r(i)*t) ;                   %solving angular displacement
    theta_prime = amplitude(i)*r(i)*cos(r(i)*t) ;            %solving angular velocity
    
    
    %graph plot
    
    
    plot(r, amplitude) ;
    
end

Unrecognized function or variable 'odefcn_script'.

Error in solution (line 19)
[t,y] = ode45(@(t,y) odefcn_script(t,y,f,r(i),zeta), tspan, Y0);

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 106)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

ODE SCRIPT

function dydt = odefcn(t,y,f,r,zeta)         % defining the equation
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = f.*sin(r.*t)-2.*zeta.*y(2)-y(1); 
end

[Steven Lord]

What specifically does "is not working" mean in this context?

• Do you receive warning and/or error messages? If so the full and exact text of those messages (all the text displayed in orange and/or red in the Command Window) may be useful in determining what's going on and how to avoid the warning and/or error.
• Does it do something different than what you expected? If so, what did it do and what did you expect it to do?
• Did MATLAB crash? If so please send the crash log file (with a description of what you were running or doing in MATLAB when the crash occured) to Technical Support using the Contact Support link on the Support section of the MathWorks website so we can investigate.

Answer 1

[Walter Roberson]

opts = odeset('RelTol',1e-6,'AbsTol',[1e-9 1e-9]);    % setting ODE solver tolerances

%% setting parameters and initial conditions

r= (0.1:0.1:1.5) ;
zeta= 0.35 ;
epsilon = -1/6 ;
Y0= [0.50, 0] ;
tspan= [0 1000] ;
f = 0.45 ;

opts = odeset('RelTol',1e-6,'AbsTol',[1e-9 1e-9]);    % setting ODE solver tolerances

%% ODE solver
for i = 1 : length(r)
    [t,y] = ode45(@(t,y) odefcn_script(t,y,f,r(i),zeta), tspan, Y0);
    
    
    harmonic_equation = @(A) f^2-(-r(i)^2*A + A +((epsilon*3*A^3)/4)^2)-((2*A*zeta*r(i))^2);
    
    amplitude(:,i) = fzero(harmonic_equation, 1);   %solving harmonic balace equation
    
    theta1 = amplitude(i)*sin(r(i)*t) ;                   %solving angular displacement
    theta_prime = amplitude(i)*r(i)*cos(r(i)*t) ;            %solving angular velocity  
end
%graph plot
    
plot(r, amplitude) ;

function dydt = odefcn_script(t,y,f,r,zeta)         % defining the equation
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = f.*sin(r.*t)-2.*zeta.*y(2)-y(1); 
end

[Alastair Strachan]

thank you so much!

Answer 2

[Cris LaPierre]

There are a few items I see.

1. You call your odefcn using odefcn_script but you named it odefcn. You have to use the same name.
2. You solve your ode for r(i), but try to plot the result against the full vector r. I suggest moving your plotting code outside your for loop.

Here's how I might modify your code

opts = odeset('RelTol',1e-6,'AbsTol',[1e-9 1e-9]);    % setting ODE solver tolerances

%% setting parameters and initial conditions
r= (0.1:0.1:1.5) ;
zeta= 0.35 ;
epsilon = -1/6 ;
Y0= [0.50, 0] ;
tspan= [0 1000] ;
f = 0.45 ;

%% ODE solver
for i = (1:1:15) ;    
    [t,y] = ode45(@(t,y) odefcn(t,y,f,r(i),zeta), tspan, Y0);
        
    harmonic_equation = @(A) f^2-(-r(i)^2*A + A +((epsilon*3*A^3)/4)^2)-((2*A*zeta*r(i))^2);
    amplitude(:,i) = fzero(harmonic_equation, 1);   %solving harmonic balace equation
    theta1 = amplitude(i)*sin(r(i)*t) ;                   %solving angular displacement
    theta_prime = amplitude(i)*r(i)*cos(r(i)*t) ;            %solving angular velocity
        
end
%graph plot
plot(r,amplitude)

function dydt = odefcn(t,y,f,r,zeta)         % defining the equation
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = f.*sin(r.*t)-2.*zeta.*y(2)-y(1); 
end