Variable is in a function file, how do I use a loop to change the value of this variable?

If the variable is in the main program, I can compile a loop to change the value of a variable and get the result with the variable, but what if the variable is in the function file?

it means i want to use loop to get the plot that zmax changes as the cd or cm changes.

main program are as this

clc;
clear all;
close all;

step=0.1;
t0=50;
length1=t0/step;

opts = odeset('RelTol',1e-4,'AbsTol',1e-4);
Z=zeros(length1,2);

[T,Z]=ode45('vdp', [0:step:t0], [0,1], opts);
 
figure(1)
plot(T,Z(:,1),'r.');

zend=length(Z(:,1));

zmax=max(Z(20:zend,1))%get the plot that zmax changes as the cd or cm changes

function file is as this

function dz = vdp(t,z)
 
g=9.8;
gamma=10094;
d=20;
J0=15;
J1=5;
c1=7680;
r=0.21;
m=13500;
G=m*g;
c2=32000;
m1=13500;
m2=0;
%************************************************************
%************************************************************
cd=0.8;%which I want to use the loop to change
%************************************************************
%************************************************************
k0=2000;
p0=1000;
area=15.9;
%************************************************************
%************************************************************
ca=1;%which I want to use the loop to change
%************************************************************
%************************************************************
h=0.5;
t0=3.5;
d1=0.8283479061205706; 
height=1;
rho=1025;
 
dz=zeros(2,1);
dz(1)=z(2);
if z(2)>=0
    a1x=1;
else a1x=0;
end
if z(2)>=0
    Jx=J0;
else Jx=J1;
end
if z(2)>=0
    b1x=0;
else b1x=0;
end
 
dz(2)=...
 ...
-(p0/3 + g*m + (k0*z(1))/9 - area*g*rho*(d1 - z(1) + (h*cos((2*pi*t)/t0))/2) + (a1x*c2*z(2))/r + (area*cd*rho*abs(z(2) + (pi*h*sin((2*pi*t)/t0))/t0)*(z(2) + (pi*h*sin((2*pi*t)/t0))/t0))/2 + (2*pi^2*area*h*rho*cos((2*pi*t)/t0)*(d1 - z(1) + (h*cos((2*pi*t)/t0))/2))/t0^2 + (2*pi^2*area*ca*h*rho*cos((2*pi*t)/t0)*(d1 - z(1) + (h*cos((2*pi*t)/t0))/2))/t0^2)/(m + Jx/r^2 + area*ca*rho*(d1 - z(1) + (h*cos((2*pi*t)/t0))/2));


end

Answer 1

[Jan]

Do not redefine the important Matlab function cd as a variable, because this can cause serious troubles during debugging.

[T,Z]=ode45('vdp', [0:step:t0], [0,1], opts);

Providing a function to be integrated as a char vector is outdated for over 15 years now. Although this is still supported for backwar compatibility, it is less powerful and secure than a function handle. In addition a function handle allows to define a parameter:

p = 3.1415;  % Your parameter
[T,Z] = ode45(@(t,y) vdp(t, y, p), 0:step:t0, [0,1], opts);

function dy = vdp(t, y, p)
...
end

See Answers: Anonymous for params

By the way, o:step:to is a vector already. Using the Matlab operator for array concatenation [] is a waste of time here only.

Another problem can cause invalid outputs: Your function to be integrated is not smooth:

if z(2)>=0
    Jx=J0;
else Jx=J1;
end

Matlab's ODE integrators are designed to process smooth functions only. Otherwise the step size controller can drive completely mad, which can cause different effects: A stop with an error message (minimal step size exceeded), or even worse a result, which is dominated by rounding errors due to a massive increase of the evaluated steps. Don't do this. Do not use a numerical tool for a computation, which does not match the specified limitations.

Use event functions instead, which stop the integration when the state of the function changes.

See ODE45 discontinuities

[dcydhb dcydhb]

why we can't use the

if z(2)>=0
    Jx=J0;
else Jx=J1;
end

and can you clarify the event functions usage by an example?

thanks a lot!!!

Answer 2

[Rik]

According to the doc you can pass extra parameters by setting an anonymous function. I see no indication in the R2014a doc that this wouldn't work in that release.

If you make the changes indicated below, you should be able to have cd and ca as inputs, so you can change them in a loop.

cd=0.7;
ca=1.1;
[T,Z]=ode45(@(t,z) vdp(t,z,cd,ca), [0:step:t0], [0,1], opts);

function dz = vdp(t,z,cd,ca)
if nargin<4
    %set defaults
    ca=1;
    if nargin<3
        cd=0.8;
    end
end

%%
%rest of your function
%%
end

[dcydhb dcydhb]

Can you be more specific?thanks！

[Rik]

See my editted answer.