How can I plot frequency response for a vibratory system with nonlinear differential equations?

Question

I have a system of nonlinear differential equations for a 3 degree of freedom vibratory system. system of differential equations

First I want to plot y, y_L and y_R against time (for a given value for Omega) and then I want to plot the domains (max values of y, y_L and y_R) against various amounts of Omega. Unfortunately, I am not good at Octave. I have written the following code in Octave (based on a sample given by one of the users), but it ends with this error: "anonymous function bodies must be single expressions".

I would be grateful if anyone can help me.

Here is the code:

Me = 4000;
me = 20;
c = 2000;
c1 = 700;
c2 = 700;
k = 20000;
k1 = 250000;
k2 = 20000;
a0 = 0.01;
om = 25;
mu1 = (c+2*c2)/(Me);
mu2 = (c2)/(Me);
mu3 = (c1+c2)/(me);
mu4 = (c2)/(me);
w12 = (2*k2)/(Me);
w22 = (k1+k2)/(me);
a1 = (k2)/(me);
a2 = (k)/(Me);
F0 = (k1*a0)/(Me);
couplode = @(t,y) [y(2); mu4*y(4) - mu3*y(2) - w22*y(1) + a1*y(3) + F0*cos(om*t); y(4); mu2*(y(2)+y(6)) - mu1*y(4) - w12*y(3) + 0.5*w12*(y(1)+y(5)) + a2((y(3)).^3; y(6); mu4*y(4) - mu3*y(6) - w22*y(5) + a1*y(3) + F0*cos(om*t)];

[t,y] = ode45(couplode, [0 0.49*pi], [1;1;1;1;1;1]*1E-8); 

figure(1)
plot(t, y)
grid
str = {'$$ \dot{y_L} $$', '$$ y_L $$', '$$ \dot{y} $$', '$$ y $$', '$$ \dot{y_R} $$', '$$ y_R $$'};
legend(str, 'Interpreter','latex', 'Location','NW')

Answer 1

You have a strange term rather at the end of the vector definition

... + a2((y(3)).^3

You certainly meant

... + a2*y(3).^3

You get better visibility and easier debugging by breaking that into separate lines

couplode = @(t,y) [ y(2); 
                    mu4*y(4)-mu3*y(2)-w22*y(1)+a1*y(3)+F0*cos(om*t); 
                    y(4); 
                    mu2*(y(2)+y(6)) - mu1*y(4) - w12*y(3) + 0.5*w12*(y(1)+y(5)) + a2*y(3).^3; 
                    y(6); 
                    mu4*y(4)-mu3*y(6)-w22*y(5)+a1*y(3)+F0*cos(om*t)];

At least in this form, spaces or no spaces makes no difference. In general in matlab/octave [a +b -c] is the same as [a, +b, -c], so one has to be careful that the expression is not interpreted as matrix row. Spaces on both sites of the operation sign switches back to the single-expression interpretation.