Problem with EventFunction triggering MATLAB

Hello! Feel free to ask me anything if it helps better understand the code. Thank you for the help.

I'm compiling this code to simulate the behaviour of a mechanical machine. The simulation is based on two odes (eq1 and eq2) and it works perfectly fine for certain values of Casp and Cman. However, when I change them to higher values, the simulation starts having issues. The simulation starts with eq1 and the code correctly changes to eq2. However, when EventFunction2 is triggered, the simulation starts going in a loop where it is continuously switching between eq1 and eq2. When I plot the results, it looks like the simulation gets stuck in eq2.

Here is the code:

%Parametri Da Cambiare
nciclo = 10; %tempo
sigma = 0.02; % 1% del volume del pistone
Casp = 30; %coefficiente aspirazione
Cman = 90; %coefficiente mandata
coeff3_out = 100; %questo va aumentato se la pressione raggiunta è troppo alta (avendo modificato sigma)
h = 0; %1E-6; %0 %spessore trafilamento %con trafilamento fa ancora fatica
ratio_volume = 0.02; %Bunsen coefficient
type = 'I'; %A 

%PARAMETRI PISTONE
N = 420; %rpm
freq = N/60; %frequenza, giri al secondo
omega = 2*pi*N/60; %rad/s
d = 0.0605; %m %diametro pistone
c = 6*d; %corsa
r = c/2; %raggio di manovella
l_pist = d; %lunghezza pistone
lambda = 1/0.1875;
Vcil = pi/4*d^2*c; %m^3 
%sfasa pistone
sfasa =  fzero(@(tett) r*(1-cos(tett)+lambda-sqrt(lambda^2-sin(tett).^2))-c/2,pi/2);

%CARTER
P_cart = 1*10^5; %1 bar

beta_0 = 1700*10^6; %Pa %rigidità olio 
rho_oil = 862; %kg/m3 %850
visc_cin = 68*10^(-6);
visc_din = rho_oil * visc_cin;
V03 = Vcil*3; %300 di olio

%PIASTRA FORATA
Nf = 30; %numero di fori
Df = 0.01; %diametro del foro
Cd = 0.8; %possibile cambiarlo 

%MEMBRANA
delta_m = 0.006; %m
r_m = sqrt(Vcil/(pi/3*2*delta_m)); %raggio
rhom = 7855.6; %kg/m3 %densità acciaio
thic = 1.5*10^-3; %m 
young_m = 210*10^9; %Pa
poisson_m = 0.3;
m_m = rhom*thic*pi*r_m^2; %massa membrana
k_m = 16*young_m/(3*(1-poisson_m^2))*thic^3*pi/((r_m)^2);
damp = 0;

%GAS
Pasp = Casp*10^6;  % Pasp = 10 MPa
Tasp = 298.15;  %K  Tamb
Pman = Cman*10^6; % Pa %pressione di scarico % Pman = 30 MPa
A_in_gas = 2*10^-4; %m^2
A_out_gas = 0.7*10^-4; %m^2 %area ragionevole
fluido = 'hydrogen';%'air.mix';

%VALVOLE OLIO
%plima = 0.01*10^5; % 0.01 bar ovvero 1 kPa
p_cp   = Pasp; 
p_over = Pman*1.1; %oil pressure should be 10% higher than gas discharge pressure
pout  = 1*10^5;

%FONDO CORSA

Unrecognized function or variable 'Vcil'.

[gamma,cmax,tetamax,deltamax] = calc_gamma(sigma,Vcil,c,d,r,lambda,r_m);
deltamax = round(deltamax,6);
delta_mezz = deltamax;
Vgas_cavity = pi/3*r_m^2*deltamax; %è con deltamax %prof
V02 = 0.2*Vgas_cavity;
V01 = 0.01*Vgas_cavity;

%Compensating Pump
c_cp = 0.01; %m 
[r_cp] = volume_CP_integrale_test(sigma,r,c_cp,Vcil,omega,c,lambda); 
%m %questo è da cambiare ogni volta
%r_cp = 0.07; %m %0.0134
A_cp = pi*r_cp^2;
delta = 0;
lambda_cp = 0.2; %NON cambia anche lambda se cambio raggio, perchè il raggio che considera in questa formula è quello di manovella

%CONDIZIONI INIZIALI (parto da membrana in quiete in mezzeria)
P3_0 = Pasp; 
P2_0 = Pasp; 
P0 = Pasp; 
T0 = Tasp; 
M0 = refpropm('D','T',T0,'P',Pasp/1000,fluido) * (V01+Vgas_cavity/2); 
U0 = M0*refpropm('U','T',Tasp,'P',Pasp/1000,fluido); %KJ %non è diviso solo per 100?? 
delta0 = 0;
deltap0 = 0;

%parameters related to odes
y0 = [P3_0 P2_0 delta0 deltap0 M0 U0];
y = y0;
tciclo = nciclo/freq;
tspan = [0 tciclo];
tstart = tspan(1);
t = tstart;
tend = tspan(end);
teout = [];     
yeout = [];    
ieout = [];
fcn = @eq1_0611_betaNC;
extdfunc1 = @(t,y)EventFunction1(t,y,deltamax);
extdfunc2 = @(t,y)EventFunction2(t,y,Pman);
opt = odeset('Events',extdfunc1);
flag = 0; %definisco la flag iniziale siccome sto iniziando con eq1
%ODES
while(t(end) < tend)
    [at,ay,ate,aye,aie] = ode23t(@(tt,y) fcn(tt,y,d,omega,lambda,r,beta_0,type,rho_oil,Cd,r_m, ...
        damp,m_m,k_m,Pasp,Tasp,Pman,A_in_gas,A_out_gas,Nf,Df, ...
        fluido,V01,V02,V03,P_cart,h,visc_din,l_pist,p_cp,p_over,pout,sfasa,deltamax,...
        coeff3_out,lambda_cp,A_cp,c_cp,delta,ratio_volume), [t(end) tend], y(end,:) , opt);
    %conserva val
    t = [t; at(2:end)];
    y = [y; ay(2:end,:)]; %mi salva i valori ogni qualvolta una delle funzioni raggiunge la sua event function
    if flag == 2 || flag == 0
        y(end,3) = deltamax;     %0.0053;
        y(end,4) = -0.00001*y(end,4); %coeff = 0.5
    end
    if flag == 1
        y(end,4) = abs(y(end,4)); 
        %impongo anche che la P2 e la P3 devono essere molto simili da quando
        %ripartono
        y(end,1) = Pman; %li ho aggiunti
        y(end,2) = Pman;
    end
    teout = [teout; ate];
    yeout = [yeout; aye];
    ieout = [ieout; aie];
    %y(end,2)
    %y(end,3)
    %y(end,4)
    %decide new equation and event function
    if y(end,4) > 0 %controllare condizioni
        fcn = @eq1_0611_betaNC;
        opt = odeset('Events', extdfunc1);
        disp('no')
        flag = 2;
    else  
        fcn = @eq2_0611_betaNC;
        opt = odeset('Events', extdfunc2);
        disp('si')
        flag = 1;
    end
end

figure(5)
plot(t,y(:,3),t,y(:,4))
legend('spostamento membrana','velocità membrana','location','southeast')

And the event functions are:

function [position,isterminal,direction] = EventFunction1(t,y,deltamax)
position = y(3)-deltamax-1e-7; % importante che non sia proprio 0.005 se no ho problemi con if
isterminal = 1;  % Halt integration 
direction = +1;   % The zero can be approached from either direction
end
function [position,isterminal,direction] = EventFunction2(t,y,Pman)
position = y(1)-Pman-10; %potrei volere che la P2 debba arrivare alle Pman
isterminal = 1;  % Halt integration 
direction = -1; %-1   % The zero can be approached from either direction
end

[Torsten]

Please prepare your code so that it can be run (see above for a missing variable) and that it incorporates the problem you are facing.

[Andrea Martella]

I added what needed. Attached are also the files needed to work the code. Let me know if you encounter any other issue.

%Parametri Da Cambiare
nciclo = 10; %tempo
sigma = 0.02; % 1% del volume del pistone
Casp = 30; %coefficiente aspirazione
Cman = 90; %coefficiente mandata
coeff3_out = 100; %questo va aumentato se la pressione raggiunta è troppo alta (avendo modificato sigma)
h = 0; %1E-6; %0 %spessore trafilamento %con trafilamento fa ancora fatica
ratio_volume = 0.02; %Bunsen coefficient
type = 'I'; %A 
%PARAMETRI PISTONE
N = 420; %rpm
freq = N/60; %frequenza, giri al secondo
omega = 2*pi*N/60; %rad/s
d = 0.0605; %m %diametro pistone
c = 6*d; %corsa
r = c/2; %raggio di manovella
l_pist = d; %lunghezza pistone
lambda = 1/0.1875;
Vcil = pi/4*d^2*c; %m^3 
%sfasa pistone
sfasa =  fzero(@(tett) r*(1-cos(tett)+lambda-sqrt(lambda^2-sin(tett).^2))-c/2,pi/2);
%CARTER
P_cart = 1*10^5; %1 bar
beta_0 = 1700*10^6; %Pa %rigidità olio 
rho_oil = 862; %kg/m3 %850
visc_cin = 68*10^(-6);
visc_din = rho_oil * visc_cin;
V03 = Vcil*3; %300 di olio
%PIASTRA FORATA
Nf = 30; %numero di fori
Df = 0.01; %diametro del foro
Cd = 0.8; %possibile cambiarlo 
%MEMBRANA
delta_m = 0.006; %m
r_m = sqrt(Vcil/(pi/3*2*delta_m)); %raggio
rhom = 7855.6; %kg/m3 %densità acciaio
thic = 1.5*10^-3; %m 
young_m = 210*10^9; %Pa
poisson_m = 0.3;
m_m = rhom*thic*pi*r_m^2; %massa membrana
k_m = 16*young_m/(3*(1-poisson_m^2))*thic^3*pi/((r_m)^2);
damp = 0;
%GAS
Pasp = Casp*10^6;  % Pasp = 10 MPa
Tasp = 298.15;  %K  Tamb
Pman = Cman*10^6; % Pa %pressione di scarico % Pman = 30 MPa
A_in_gas = 2*10^-4; %m^2
A_out_gas = 0.7*10^-4; %m^2 %area ragionevole
fluido = 'hydrogen';%'air.mix';
%VALVOLE OLIO
%plima = 0.01*10^5; % 0.01 bar ovvero 1 kPa
p_cp   = Pasp; 
p_over = Pman*1.1; %oil pressure should be 10% higher than gas discharge pressure
pout  = 1*10^5;
%FONDO CORSA
[gamma,cmax,tetamax,deltamax] = calc_gamma(sigma,Vcil,c,d,r,lambda,r_m);
deltamax = round(deltamax,6);
delta_mezz = deltamax;
Vgas_cavity = pi/3*r_m^2*deltamax; %è con deltamax %prof
V02 = 0.2*Vgas_cavity;
V01 = 0.01*Vgas_cavity;
%Compensating Pump
c_cp = 0.01; %m 
[r_cp] = volume_CP_integrale_test(sigma,r,c_cp,Vcil,omega,c,lambda); 
%m %questo è da cambiare ogni volta
%r_cp = 0.07; %m %0.0134
A_cp = pi*r_cp^2;
delta = 0;
lambda_cp = 0.2; %NON cambia anche lambda se cambio raggio, perchè il raggio che considera in questa formula è quello di manovella
%CONDIZIONI INIZIALI (parto da membrana in quiete in mezzeria)
P3_0 = Pasp; 
P2_0 = Pasp; 
P0 = Pasp; 
T0 = Tasp; 
M0 = refpropm('D','T',T0,'P',Pasp/1000,fluido) * (V01+Vgas_cavity/2); 

Unrecognized function or variable 'refpropm'.

U0 = M0*refpropm('U','T',Tasp,'P',Pasp/1000,fluido); %KJ %non è diviso solo per 100?? 
delta0 = 0;
deltap0 = 0;
%parameters related to odes
y0 = [P3_0 P2_0 delta0 deltap0 M0 U0];
y = y0;
tciclo = nciclo/freq;
tspan = [0 tciclo];
tstart = tspan(1);
t = tstart;
tend = tspan(end);
teout = [];     
yeout = [];    
ieout = [];
fcn = @eq1_0611_betaNC;
extdfunc1 = @(t,y)EventFunction1(t,y,deltamax);
extdfunc2 = @(t,y)EventFunction2(t,y,Pman);
opt = odeset('Events',extdfunc1);
flag = 0; %definisco la flag iniziale siccome sto iniziando con eq1
%ODES
while(t(end) < tend)
    [at,ay,ate,aye,aie] = ode23t(@(tt,y) fcn(tt,y,d,omega,lambda,r,beta_0,type,rho_oil,Cd,r_m, ...
        damp,m_m,k_m,Pasp,Tasp,Pman,A_in_gas,A_out_gas,Nf,Df, ...
        fluido,V01,V02,V03,P_cart,h,visc_din,l_pist,p_cp,p_over,pout,sfasa,deltamax,...
        coeff3_out,lambda_cp,A_cp,c_cp,delta,ratio_volume), [t(end) tend], y(end,:) , opt);
    %conserva val
    t = [t; at(2:end)];
    y = [y; ay(2:end,:)]; %mi salva i valori ogni qualvolta una delle funzioni raggiunge la sua event function
    if flag == 2 || flag == 0
        y(end,3) = deltamax;     %0.0053;
        y(end,4) = -0.00001*y(end,4); %coeff = 0.5
    end
    if flag == 1
        y(end,4) = abs(y(end,4)); 
        %impongo anche che la P2 e la P3 devono essere molto simili da quando
        %ripartono
        y(end,1) = Pman; %li ho aggiunti
        y(end,2) = Pman;
    end
    teout = [teout; ate];
    yeout = [yeout; aye];
    ieout = [ieout; aie];
    %y(end,2)
    %y(end,3)
    %y(end,4)
    %decide new equation and event function
    if y(end,4) > 0 %controllare condizioni
        fcn = @eq1_0611_betaNC;
        opt = odeset('Events', extdfunc1);
        disp('no')
        flag = 2;
    else  
        fcn = @eq2_0611_betaNC;
        opt = odeset('Events', extdfunc2);
        disp('si')
        flag = 1;
    end
end
figure(5)
plot(t,y(:,3),t,y(:,4))
legend('spostamento membrana','velocità membrana','location','southeast')

function [position,isterminal,direction] = EventFunction1(t,y,deltamax)
position = y(3)-deltamax-1e-7; % importante che non sia proprio 0.005 se no ho problemi con if
isterminal = 1;  % Halt integration 
direction = +1;   % The zero can be approached from either direction
end
function [position,isterminal,direction] = EventFunction2(t,y,Pman)
position = y(1)-Pman-10; %potrei volere che la P2 debba arrivare alle Pman
isterminal = 1;  % Halt integration 
direction = -1; %-1   % The zero can be approached from either direction
end
function [] = test()
%clc 
%clear
 
%d =  0.0605; %diametro del pistone %m 
%c = 6*d; %corsa 
%V = pi/4*d^2*c; %m^3
%r = c/2; %raggio di manovella
%lambda = 1/0.1875; %4 %rapporto tra connectin rod e crank radius %scelto 

%CP
A_cp = pi*0.0182^2; %m^2
V_cp = sigma*Vcil; % m^3
lambda_cp = 0.2;
%N = 420; % rpm 
%omega = 2*pi*N/60;
%f = N/60;
%t_ciclo = 1/f; %devo moltiplicarlo per il tempo giusto di apertura che non è 0.14
%dividendo così i due contributi riesco ad avere sempre portate entranti
%positive
t_cp_1 = 0.035; %tempo per cui è aperta CP dopo la mezzeria %0.038
t_cp_2 = 0.031; %tempo per cui è aperta CP prima della mezzeria %0.018 
t_cp_tot = t_cp_1 + t_cp_2; %tempo di aspirazione calcolato dai grafici della simulazione

%per fare un conto più corretto posso trovare l'angolo a cui CP subentra
sfasa =  fzero(@(tett) r*(1-cos(tett)+lambda-sqrt(lambda^2-sin(tett).^2))-c/2,pi/2);
sfasa_cp_tot = sfasa - omega*t_cp_2; %influenza questo %io avevo messo 2*pi MA è SBAGLIATO

f = @(t,A) omega*c_cp/2*A*(sin(omega*t + sfasa_cp_tot)+lambda_cp*sin(omega*t + sfasa_cp_tot).*cos(omega*t + sfasa_cp_tot)/(sqrt(1-lambda_cp^2*sin(omega*t+sfasa_cp_tot).^2)));

%per calcolare nuova area

A_cpnew = fzero(@(A) V_cp - integral(@(t) f(t,A), 0, t_cp_tot), A_cp) ; %mettere un intervallo anzichè A_cp
r_cpnew = sqrt(A_cpnew/pi);

%funcheck = @(t) c_cp/2*A_cpnew*(sin(omega*t + sfasa_cp_1)+lambda_cp*sin(omega*t + sfasa_cp_1).*cos(omega*t + sfasa_cp_1)/(sqrt(1-lambda_cp^2*sin(omega*t+sfasa_cp_1).^2)));
%V_cp_check = integral(funcheck,0,t_cp);

fcheck = @(t) omega*c_cp/2*A_cpnew*(sin(omega*t + sfasa_cp_1)+lambda_cp*sin(omega*t + sfasa_cp_1).*cos(omega*t + sfasa_cp_1)/(sqrt(1-lambda_cp^2*sin(omega*t+sfasa_cp_1).^2)));
%gcheck = @(t) omega*c_cp/2*A_cpnew*(sin(omega*t + sfasa_cp_2)+lambda_cp*sin(omega*t + sfasa_cp_2).*cos(omega*t + sfasa_cp_2)/(sqrt(1-lambda_cp^2*sin(omega*t+sfasa_cp_2).^2)));

%V_cp_check = integral(fcheck,0,t_cp_1) + integral(gcheck,0,t_cp_2); %cos' in teoria le mie portate sono sempre positive

%t = linspace(0,t_cp_tot,100);
%figure(1)
%plot(t,fcheck(t),t,gcheck(t))


%c0 = fzero(@(c0) (1-alpha) - integral(@(y) fun(y,c0), 0,(2.*(aLehat./ c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  ))),   0.01);
end
function [gamma,cmax,tetamax,deltamax] = calc_gamma(sigma,Vcil,c,d,r,lambda,r_m)
%calcolo dei vari parametri a seconda della volume che voglio buttare fuori
%

%clc
%clear
%V_over = 1.04354*10^(-5);       %1.73922594468318*10^(-4) t.c. gamma = 0.2
%N = 420; % rpm %420
%Hz = N/60; % frequenza, giri al secondo
%omega = 2*pi*N/60; % rad/s
%d =  0.0605; %diametro del pistone %m
%c = 6*d; %corsa
%l_pist = d;
%r = c/2; %raggio di manovella
%lambda = 1/0.1875;
%r_m = 0.288;
V_over = sigma*Vcil;

cmax = fzero(@(cmax) V_over - pi*(d/2)^2*(c-cmax),c);

tetamax = fzero(@(tett) r*(1-cos(tett)+lambda-sqrt(max(0, lambda^2-sin(tett).^2))) - cmax, 2*pi);

deltamax = (cmax*pi*(d/2)^2)/(pi/3*r_m^2*2);

gamma = (c*pi*(d/2)^2)/(pi/3*r_m^2*2*deltamax)-1;

end
function dy = eq2_0611_betaNC(tt,y,d,omega,lambda,r,beta_0,type,rho_oil,Cd,r_m, ...
        damp,m_m,k_m,Pasp,Tasp,Pman,A_in_gas,A_out_gas,Nf,Df, ...
        fluido,V01,V02,V03,P_cart,h,visc_din,l_pist,p_cp,p_over,pout,sfasa,deltamax,...
        coeff3_out,lambda_cp,A_cp,c_cp,delta,ratio_volume)

teta = omega*tt+sfasa;

beta3 = beta_Cho(ratio_volume, beta_0, pout, y(1),Tasp,type); 
beta2 = beta_Cho(ratio_volume, beta_0, pout, y(2),Tasp,type); 

 x  = r*(1-cos(teta)+lambda-sqrt(lambda^2-sin(teta)^2)); %ricordati il punto nell'elevazione a potenza %
 xp = r*( sin(teta)  +0.5/sqrt(lambda^2-sin(teta)^2).*sin(2*teta))*omega; %così è rispetto il tempo siccome c'è omega

 V3    = V03 + pi/4*d^2*x; %in questo modo, il volume V3 sta aumentando, quindi il pistone sta andando verso il basso
 dV3   = pi/4*d^2*xp;

 if p_cp > y(1)
     Q3_in = omega*c_cp/2*A_cp*(sin(teta-delta)+lambda_cp*sin(teta-delta)*cos(teta-delta)/sqrt(1-lambda_cp^2*sin(teta-delta)^2));
 else
     Q3_in = 0;
 end
 if y(1) > p_over + pout
    Q3_out = - coeff3_out*(0.001/60/10^5)*(y(1)-p_over-pout);
 else
    Q3_out = 0;
 end

 Q3_ext = Q3_in + Q3_out;

  %trascinamento
 if dV3 > 0 %il pistone sta scendendo
     Q3_tras = - xp/2 * pi*h*l_pist;
 elseif dV3 == 0
     Q3_tras = 0;
 elseif dV3 < 0 
     Q3_tras =  - xp/2 * pi*h*l_pist; %ci devo mettere - siccome xp è negativo
 end
  
 %trafilamento
    Q3_traf = - (y(1)-P_cart)*pi*d*h^3/(12*visc_din*l_pist);  %esce sempre

 if y(1) > y(2)
     Q23 = - Cd * Nf * pi/4 * Df^2 * sqrt(2*(y(1)-y(2))/rho_oil); 
 else
     Q23 = Cd * Nf * pi/4 * Df^2 * sqrt(2*(y(2)-y(1))/rho_oil);
 end

 Q3 = Q23 + Q3_ext + Q3_traf + Q3_tras;

 %Area 2
 Q2 = -Q23;

 V2 = V02 + pi/3*r_m^2*(deltamax + y(3));
 dV2 = pi/3*r_m^2*y(4);

 V1 = V01 + pi/3*r_m^2*(deltamax - y(3));
 dV1 = -dV2;

 rhogas = y(5)/V1; %KG/m^3 
 ugas  = y(6)/y(5); %KJ/KG
 Pgas = refpropm('P','D',rhogas,'U',ugas,fluido)*1000; % Pa 
 
 hi  = refpropm('H','T',Tasp,'P',Pasp/1000,fluido); %sicuro che sia /1000?
 ho  = refpropm('H','D',rhogas,'U',ugas,fluido); 

  if Pasp - Pgas > 0 %processo di aspirazione
      PR_asp = Pgas/Pasp; %ATTENZIONE A PASP
      k = refpropm('K','T',Tasp,'P',Pasp/1000,fluido); 
      %perchè qui k assume un unico valore, mentre in discharge k varia con il variare della densità? --> prob perchè il discharge dipende dalle condizioni interne alla 
      %camera, mentre in aspirazione dipendono solo da condizioni esterne

      m_gin = (A_in_gas)*sqrt(Pasp*refpropm('D','T',Tasp,'P',Pasp/1000,fluido)) ... %rhogas iniziale
                        *sqrt(2/(k-1)*(PR_asp^(2/k)-PR_asp^((k+1)/k)));
 else
      m_gin = 0;
 end

 if Pgas - Pman > 0  %processo di scarico
     k  = refpropm('K','D',rhogas,'U',ugas,fluido); 
     PR_man = Pman/Pgas;

     if PR_man > (2/(k+1))^(k/(k-1)) %cioè del rapporto critico di pressione %è divergente o convergente questo ugello??
          m_gout = (A_out_gas) *sqrt(Pgas*rhogas) ...
                        *sqrt(2/(k-1)*(PR_man^(2/k)-PR_man^((k+1)/k)));
     else
          m_gout = (A_out_gas) *sqrt(Pgas*rhogas) ...
                        *sqrt(k*(2/(k+1))^((k+1)/(k-1)));  
                                                        %formula usata nel caso in cui l'ugello convergente-divergente è critico --> la portata non dipende più dal PR                                                   
     end    
 else
          m_gout = 0; 
 end

 dy(1) = (beta3/V3)*(Q3-dV3); 
 dy(2) = (beta2/V2)*(Q2-dV2); 
 dy(3) = 0;
 dy(4) = 0; %(y(2)-Pgas)/m_m*pi*r_m^2; %devo mettere il valore assoluto? %prof
 dy(5) = (m_gin-m_gout);
 dy(6) = (m_gin*hi - m_gout*ho) - Pgas*dV1; %formula adiabatica (e senza perdite di olio)

 dy = dy(:);
end

function dy = eq1_0611_betaNC(tt,y,d,omega,lambda,r,beta_0,type,rho_oil,Cd,r_m, ...
        damp,m_m,k_m,Pasp,Tasp,Pman,A_in_gas,A_out_gas,Nf,Df, ...
        fluido,V01,V02,V03,P_cart,h,visc_din,l_pist,p_cp,p_over,pout,sfasa,deltamax,...
        coeff3_out,lambda_cp,A_cp,c_cp,delta,ratio_volume)

teta = omega*tt+sfasa;

beta3 = beta_Cho(ratio_volume, beta_0, pout, y(1),Tasp,type);
beta2 = beta_Cho(ratio_volume, beta_0, pout, y(2),Tasp,type); 

 x  = r*(1-cos(teta)+lambda-sqrt(lambda^2-sin(teta)^2)); %ricordati il punto nell'elevazione a potenza %
 xp = r*( sin(teta)  +0.5/sqrt(lambda^2-sin(teta)^2).*sin(2*teta))*omega; %così è rispetto il tempo siccome c'è omega

 V3    = V03 + pi/4*d^2*x; %in questo modo, il volume V3 sta aumentando, quindi il pistone sta andando verso il basso
 dV3   = pi/4*d^2*xp;

 if p_cp > y(1)
     Q3_in = omega*c_cp/2*A_cp*(sin(teta-delta)+lambda_cp*sin(teta-delta)*cos(teta-delta)/sqrt(1-lambda_cp^2*sin(teta-delta)^2));
 else
     Q3_in = 0;
 end
 if y(1) > p_over + pout
    Q3_out = - coeff3_out*(0.001/60/10^5)*(y(1)-p_over-pout);
 else
    Q3_out = 0;
 end

 Q3_ext = Q3_in + Q3_out;

 %trascinamento
 if dV3 > 0 %il pistone sta scendendo
     Q3_tras = - xp/2 * pi*h*l_pist;
 elseif dV3 == 0
     Q3_tras = 0;
 elseif dV3 < 0 
     Q3_tras =  - xp/2 * pi*h*l_pist; %ci devo mettere - siccome xp è negativo
 end
  
 %trafilamento
    Q3_traf = - (y(1)-P_cart)*pi*d*h^3/(12*visc_din*l_pist);  %esce sempre

 if y(1) > y(2)
     Q23 = - Cd * Nf * pi/4 * Df^2 * sqrt(2*(y(1)-y(2))/rho_oil); 
 else
     Q23 = Cd * Nf * pi/4 * Df^2 * sqrt(2*(y(2)-y(1))/rho_oil);
 end

 Q3 = Q23 + Q3_ext + Q3_tras + Q3_traf;

 %Area 2
 Q2 = -Q23;

 V2 = V02 + pi/3*r_m^2*(deltamax + y(3));
 dV2 = pi/3*r_m^2*y(4);

 V1 = V01 + pi/3*r_m^2*(deltamax - y(3));
 dV1 = -dV2;

 rhogas = y(5)/V1; %KG/m^3 
 ugas  = y(6)/y(5); %KJ/KG
 Pgas = refpropm('P','D',rhogas,'U',ugas,fluido)*1000; % Pa 
 
 hi  = refpropm('H','T',Tasp,'P',Pasp/1000,fluido); %sicuro che sia /1000?
 ho  = refpropm('H','D',rhogas,'U',ugas,fluido); 

  if Pasp - Pgas > 0 %processo di aspirazione
      PR_asp = Pgas/Pasp; %ATTENZIONE A PASP
      k = refpropm('K','T',Tasp,'P',Pasp/1000,fluido); 
      %perchè qui k assume un unico valore, mentre in discharge k varia con il variare della densità? --> prob perchè il discharge dipende dalle condizioni interne alla 
      %camera, mentre in aspirazione dipendono solo da condizioni esterne

      m_gin = (A_in_gas)*sqrt(Pasp*refpropm('D','T',Tasp,'P',Pasp/1000,fluido)) ... %rhogas iniziale
                        *sqrt(2/(k-1)*(PR_asp^(2/k)-PR_asp^((k+1)/k)));
 else
      m_gin = 0;
 end

 if Pgas - Pman > 0  %processo di scarico
     k  = refpropm('K','D',rhogas,'U',ugas,fluido); 
     PR_man = Pman/Pgas;

     if PR_man > (2/(k+1))^(k/(k-1)) %cioè del rapporto critico di pressione %è divergente o convergente questo ugello??
          m_gout = (A_out_gas) *sqrt(Pgas*rhogas) ...
                        *sqrt(2/(k-1)*(PR_man^(2/k)-PR_man^((k+1)/k)));
     else
          m_gout = (A_out_gas) *sqrt(Pgas*rhogas) ...
                        *sqrt(k*(2/(k+1))^((k+1)/(k-1)));  
                                                        %formula usata nel caso in cui l'ugello convergente-divergente è critico --> la portata non dipende più dal PR                                                   
     end    
 else
          m_gout = 0; 
 end

 dy(1) = (beta3/V3)*(Q3-dV3); 
 dy(2) = (beta2/V2)*(Q2-dV2); 
 dy(3) = y(4);
 dy(4) = -(damp/m_m)*y(4) - (k_m/m_m)*y(3) + (y(2)-Pgas)/m_m*pi*r_m^2; %devo mettere il valore assoluto? %prof
 dy(5) = (m_gin-m_gout);
 dy(6) = (m_gin*hi - m_gout*ho) - Pgas*dV1; %formula adiabatica (e senza perdite di olio)
 
 %y(1)
 dy = dy(:);
end
function beta = beta_Cho(ratio_volume,beta_0,pout,P,Tasp,type)
%clc
%clear

%Pmax = 100*10^6;
%Tasp = 300; %k
%beta_0 = 1300*10^6;
%ratio_volume = 0.09; %VG0/VLO
%pout = 101325;
if type == 'I'
    k = 1;
end
if type == 'A'
    k = 1.4; %refpropm('K','T',Tasp,'P',P/1000,'air.mix'); %qual è l'altra variabile?
end
beta = (beta_0).*(ratio_volume + (P./pout).^(1/k))./((P./pout).^(1/k)+ ratio_volume/k * beta_0./P);

%fun = @(P) (beta_0).*(ratio_volume + (P./pout).^(1/k))./((P./pout).^(1/k)+ ratio_volume/k * beta_0./P);
%P = linspace(0, Pmax, 1000);
%y = fun(P);

%figure(1)
%plot(P, y)
end

[Torsten]

Please hardcode the values coming from refprop - refpropm is not known.

[Andrea Martella]

I am sorry, but I'm not able to do that. I can attach the images of the results I obtain. Otherwise, I don't know how to help you

[Torsten]

I doubt that anybody will then be able to help you because we cannot execute your code.

[Andrea Martella]

Alright, thank you anyways. I could provide you the vectors of the parameters calculated by refprop, but more than that I really don't know how to help.