How can i plot the response of function as a vector ?

Hi,

To resolve a secon order differential equation i do this function , but i have a problem when i want to plot xdot.

function xdot = equacte(t,x)
% Function file for mass with spring.
% Position is first variable, velocity is second variable,load is the third variable
freq=100; %frequency (Hz)
w=2*pi*freq;
m=0.0112;
k=262.8257 ;
teta =-4.40e+03;
cp= 1.8492*10^-7;
for rload =0:20:50000
   A = [0,1,0;-k/m,0,-teta/m;(-teta/rload),0,-1/(rload*cp)];
   B = [0;1/m;0];
   f =sin(w*t);
   xdot = A*x+B*f;
end
end

In the principle program :

[t,x]=ode45(@equacte,[0:0.0004:1],[0,0,0]);
xdot=equacte(t,x);
plot(t,xdot(:,1),'O',t,xdot(:,2),'r',t,xdot(:,3),'b');

Answer 1

[Azzi Abdelmalek]

What this line is doing? xdot=equacte(t,x) just write

[t,x]=ode45(@equacte,[0:0.0004:1],[0,0,0]);
plot(t,x(:,1),'O',t,x(:,2),'r',t,x(:,3),'b')

[Mallouli Marwa]

But, i want to plot xdot not x

[Azzi Abdelmalek]

Now I noticed, the for loop inside your function is useless,

for rload =0:20:50000
 A = [0,1,0;-k/m,0,-teta/m;(-teta/rload),0,-1/(rload*cp)];
 B = [0;1/m;0];
 f =sin(w*t);
 xdot = A*x+B*f;
end

will produce the same effect as

   rload =50000
   A = [0,1,0;-k/m,0,-teta/m;(-teta/rload),0,-1/(rload*cp)];
   B = [0;1/m;0];
   f =sin(w*t);
   xdot = A*x+B*f;

[Azzi Abdelmalek]

If you want to plot xdot:

xdot=diff(x)
plot(t,xdot(:,1),'O',t,xdot(:,2),'r',t,xdot(:,3),'b')

[Mallouli Marwa]

No, tis for loop is necessary. For each vlue of rload i want the value of x.

[Azzi Abdelmalek]

Ok, but like you did it, the loop has no effect, only the last value will be taken in account.

[Azzi Abdelmalek]

Your function should look like

function xdot = equacte(t,x,rload)
% Function file for mass with spring.
% Position is first variable, velocity is second variable,load is the third variable
freq=100; %frequency (Hz)
w=2*pi*freq;
m=0.0112;
k=262.8257 ;
teta =-4.40e+03;
cp= 1.8492*10^-7;
A = [0,1,0;-k/m,0,-teta/m;(-teta/rload),0,-1/(rload*cp)]
B = [0;1/m;0];
f =sin(w*t);
xdot = A*x+B*f;
end

Then call it in the loop

for rload=2000:20:2040
  [t,x]=ode45(@(t,x) equacte(t,x,rload),[0 1],[0,0,0]);
  figure;
  plot(t(2:end-1),diff(x))
end