how to extract an extra variable which is calculated in ode function to main workspace?

I have to extract the values of "p" variable which is calculated in 'model' function. can anyone help?

 clear;
clc;
ifp = fopen('loma.txt','r');
 
% capture and discard the first two lines of the file
for i=1:46
   
R = fgetl(ifp); % delete redundant lines
end
 
 
% load the acceleration data 90 degree
data90 = fscanf(ifp,'%f')/100; % m/s/s
% close the input file
fclose(ifp);
ifp = fopen('loma.txt','r');
 
% capture and discard the first two lines of the file
for i=1:1646
   
F = fgetl(ifp); % delete redundant lines
end

 
% load the acceleration data 0 degree
data0 = fscanf(ifp,'%f')/100; % m/s/s
% close the input file
fclose(ifp);

dt = 0.02;
npoint=length(data90);
time=0:dt:(npoint-1)*dt;
t=time';
kx=11*10^4; %N/m
ky=8*10^4;  %N/m
fyxp=550*10^0; %N
fyxn=650*10^0; %N
fyyp=440*10^0; %N
fyyn=520*10^0; %N
alphax=0.01;
alphay=0.01;
zux0p=fyxp/kx;
zuy0p=fyyp/ky;
zux0n=fyxn/kx;
zuy0n=fyyn/ky;

Z0(1)=0;
Z0(2)=0;
Z0(3)=0;
Z0(4)=0;
Z0(5)=0;
Z0(6)=0;
[t,Z]=ode23s(@model,t,Z0,[],data0,data90,t);
qx=alphax*kx*Z(:,1)+(1-alphax)*ky*Z(:,5);
qy=alphay*ky*Z(:,3)+(1-alphay)*zuy0p/zux0p*ky*Z(:,6);


function [ydot]=model(t,Y,ax2,ay2,t1)

ax = interp1(t1,ax2,t);
ay = interp1(t1,ay2,t);
kx=11*10^4; %N/m
ky=8*10^4;  %N/m
fyxp=550*10^0; %N
fyxn=650*10^0; %N
fyyp=440*10^0; %N
fyyn=520*10^0; %N
alphax=0.01;
alphay=0.01;
beta=0.7;
dt = 0.02;
eta=1;
nu=1;
m=50;
xi=5/100;
s=8;
sigma=.2;
mu=0.1;
zux0p=fyxp/kx;
zuy0p=fyyp/ky;
zux0n=fyxn/kx;
zuy0n=fyyn/ky;
zu=((1+sign(Y(5)))*abs(zux0p)+(1-sign(Y(5)))*abs(zux0n))/2/nu;
Cx=Y(5)/zu*(1+beta*(sign(Y(2)*Y(5))-1));
Cy=Y(6)/zu*(1+beta*(sign(Y(4)*Y(6))-1));
px=sqrt(2/pi)*s/sigma*exp(-0.5*((sign(Y(2))*Y(5)/zu-mu)/sigma)^2);
py=sqrt(2/pi)*s/sigma*exp(-0.5*((sign(Y(4))*Y(6)/zu-mu)/sigma)^2);
H11=eta+(1-Y(5)/zu*Cx)*px;
H12=-Y(5)/zu*Cy*py;
H21=-Y(6)/zu*Cx*px;
H22=eta+(1-Y(6)/zu*Cy)*py;
G1=(1-Y(5)/zu*Cx)*Y(2)-Y(5)*zux0p/zu/zuy0p*Cy*Y(4);
G2=-Y(6)/zu*Cx*Y(2)+(1-Y(6)/zu*Cy)*zux0p/zuy0p*Y(4);
Ydot(5)=1/(H11*H22-H12*H21)*(H22*G1-H12*G2);
zdashx=1/(H11*H22-H12*H21)*(H22*G1-H12*G2);
Ydot(6)=1/(H11*H22-H12*H21)*(H11*G2-H21*G1);
zdashy=1/(H11*H22-H12*H21)*(H11*G2-H21*G1);
Ydot(1)=Y(2);
Ydot(3)=Y(4);
qx=alphax*kx*Y(1)+(1-alphax)*ky*Y(5);
qy=alphay*ky*Y(3)+(1-alphay)*zuy0p/zux0p*ky*Y(6);
Ydot(2)=-ax-2*xi*sqrt(kx/m)*Y(2)-qx/m;
Ydot(4)=-ay-2*xi*sqrt(ky/m)*Y(2)-qy/m;

p=1/zu^2*(Y(5)*(Y(2)-mu*zdashx)+Y(6).*(zux0p/zuy0p*Y(4)-mu*zdashy))

% v=cumtrapz(dt,p);
% Energy=(1-alphax)/2.*v;
ydot=[Ydot(1);Ydot(2);Ydot(3);Ydot(4);Ydot(5);Ydot(6)];

[Stephen23]

Just return it as an output argument, as shown here:

https://www.mathworks.com/matlabcentral/answers/512891-passing-extra-output-arguments-from-ode

https://www.mathworks.com/matlabcentral/answers/429179-how-do-i-extract-an-intermediate-variable-calculated-and-used-inside-my-ode45-function

Avoid obsolete, undocumented syntaxes by parameterizing the function properly:

https://www.mathworks.com/matlabcentral/answers/517469-some-quetions-about-ode45#answer_426013

https://www.mathworks.com/matlabcentral/answers/469314-ode45-varargin-additional-arguments-after-options#answer_381161

The recommended approach (e.g. using an anonymous function):

https://www.mathworks.com/help/matlab/math/parameterizing-functions.html

Answer 1

[Stephen23]

Ignore advice that you should use assignin or evalin, because this does not take into account how ODE solvers work, as well as being slow and complex (this topic has been discussed many times on this forum).

Here is a simple way to retrieve the p value from your function for exactly the t and Z values returned by the ODE solver.

Unlike the incomplete advice to use assignin/evalin here is a complete, simple, efficient, working solution using exactly the method that I gave you in my earlier comment. I also changed your obsolete syntax for providing extra inputs to the function with the recommended function parameterization:

https://www.mathworks.com/help/matlab/math/parameterizing-functions.html

I marked the altered lines with !!! so you can find them easily:

ifp = fopen('loma.txt','r');
% capture and discard the first two lines of the file
for i=1:46
    fgetl(ifp); % delete redundant lines
end
% load the acceleration data 90 degree
data90 = fscanf(ifp,'%f')/100; % m/s/s
% close the input file
fclose(ifp);
ifp = fopen('loma.txt','r');
% capture and discard the first two lines of the file
for i=1:1646
   fgetl(ifp); % delete redundant lines
end
% load the acceleration data 0 degree
data0 = fscanf(ifp,'%f')/100; % m/s/s
% close the input file
fclose(ifp);

dt = 0.02;
npoint=length(data90);
time=0:dt:(npoint-1)*dt;
kx=11*10^4; %N/m
ky=8*10^4;  %N/m
fyxp=550*10^0; %N
fyxn=650*10^0; %N
fyyp=440*10^0; %N
fyyn=520*10^0; %N
alphax=0.01;
alphay=0.01;
zux0p=fyxp/kx;
zuy0p=fyyp/ky;
zux0n=fyxn/kx;
zuy0n=fyyn/ky;


Z0 = [0;0;0;0;0;0];
fun = @(T,X)model(T,X,data0,data90,time(:));% !!! replace obsolete syntax !!!
[t_out,Z_out] = ode23s(fun,time(:),Z0)% !!! and use different variable names !!!

t_out = 2000×1

0 0.0200 0.0400 0.0600 0.0800 0.1000 0.1200 0.1400 0.1600 0.1800

Z_out = 2000×6

0 0 0 0 0 0 -0.0000 -0.0010 0.0000 0.0001 -0.0000 0.0000 -0.0000 -0.0010 -0.0000 -0.0005 -0.0000 -0.0000 -0.0000 0.0002 -0.0000 -0.0011 -0.0000 -0.0000 -0.0000 0.0022 -0.0000 -0.0015 -0.0000 -0.0000 0.0000 0.0044 -0.0001 -0.0019 0.0000 -0.0000 0.0002 0.0061 -0.0001 -0.0030 0.0000 -0.0000 0.0003 0.0065 -0.0002 -0.0054 0.0000 -0.0000 0.0004 0.0054 -0.0003 -0.0083 0.0000 -0.0000 0.0005 0.0032 -0.0005 -0.0104 0.0000 -0.0000

[~,p] = cellfun(fun, num2cell(t_out), num2cell(Z_out,2), 'uni',0); % !!! p values
p = cell2mat(p) % !!! and here is your output p vector !!!

p = 2000×1

0 0.0000 0.0000 0.0000 0.0000 0.0004 0.0017 0.0038 0.0061 0.0083

[subbu]

Thanks for the help. I really appreciate it

Answer 2

[Sai Veeramachaneni]

Hi,

You can leverage assignin and evalin functions for your usecase.

Hope it helps.