le fichier est le suivant:
function f=eqqdd(t,y)
global tq ca cv R r Ts j Tp alphah alphas alphap bethah bethas bethap lamdah P0
bethah=1.17;
bethas=7;
bethap=7;
alphah=0.84;
alphas=93;
alphap=93;
ca=1.55;
cv=519;
R=1.05;
V=67.5;
r=0.068;
P0=93;
j(Ts,Tp)=(alphas*Ts/1+(lamdah*Tp))-bethah*Tp;
Ts=1/(1+(tq/alphas).^bethas);
Tp=1/(1+(alphap/y1(t)).^bethap)
f = zeros(size(y),1);
f(1)=-(1/(ca*R))*y1(t)+(1/(ca*R))*y2(t)+(V/ca)*y3(t);
f(2)=(1/(cv*R))*y1(t)-((1/(cv*R))+(1/(cv*r)))*y2(t);
f(3)=j(Ts,Tp)