Would like to graph an ode45 answer with if statment

function diffy2=Ljallo_f_q2(~,X)
y=X(1);T=X(2);
%%data
Q=0.06555;           %%m^3/s
ctot=0.0203;         %%knol.m^3
alpha=26900;         %%m^2/m^3
V=6*10^-4;           %%m^3
Mhat=30;             %%kg/kmol
cpg=1070;            %%J/kg*k
delhrxn=-2.84*10^8;  %%j/kmol
ep=0.68;             %%unitless
cps=1000;            %%j/kg*k
rhog=ctot*Mhat;      %%kg/m^3
rhos=2500;           %%kg/m^3

if t<1
    yin=.01;Tin=600;
    
if (1<t)&&(t<2)
    yin=.03;
    Tin=700;
    
else
    yin=.03;
    Tin=700;
    
end

end 
zspan=[0 2];
ic=[yin Tin]
k1=6.70*10^10*exp(-12556/T);
K1=65.5*exp(961/T);
r=(0.05*k1*y)/(T(1+K1*y)^2);
diffy2(1)=(Q*ctot*(yin-y)-alpha*V*r)/(ep*V*ctot);
diffy2(2)=Q*ctot*Mhat*cpg*(Tin-T)+alpha*V*(-delhrxn)*r/(V*(ep*rhog*cpg+(1-ep)*rhos*cps));
diffy2=diffy2';

end
>> [t,y]=ode45(@Ljallo_f_q2,zspan,ic)
Undefined function or variable 't'.

Error in Ljallo_f_q2 (line 16)
if t<1

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets
args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode,
  tspan, y0, options, varargin);
 
>> 

I would like to graph this ODE with the if statemetn. The graph goes from 0 to 2 on the x axis(t) with the if conditions between 0 and 2. The initial conditions change based on the x-axis(t).I understand that my t is not defined but I dont know how to arrange it so that matlab reads it properly.

Answer 1

[Walter Roberson]

Change

function diffy2=Ljallo_f_q2(~,X)

to

function diffy2=Ljallo_f_q2(t,X)

Also you will need to move

zspan=[0 2];
ic=[yin Tin];

into your script that runs ode45 instead of having them inside your Ljallo_f_q2

Do not be surprised if ode45 gives you completely wrong answers, or complains about not being able to find a solution at time 1. The mathematics used by the ode*() functions require that your equations are continuous to two derivatives, but your equations are not continuous at t==1 or t==2.

What you should be doing is running ode45 from time 0 to 1*(1-eps), then stopping it and using the output as the initial conditions to running ode45 from time 1 to 2*(1-eps), and then stopping it and using the output as the initial conditions to run time 2 exactly.

[Rodrigo Blas]

function diffy2=Ljallo_f_q2(t,X)
y=X(1);T=X(2);
%%data
Q=0.06555;           %%m^3/s
ctot=0.0203;         %%knol.m^3
alpha=26900;         %%m^2/m^3
V=6*10^-4;           %%m^3
Mhat=30;             %%kg/kmol
cpg=1070;            %%J/kg*k
delhrxn=-2.84*10^8;  %%j/kmol
ep=0.68;             %%unitless
cps=1000;            %%j/kg*k
rhog=ctot*Mhat;      %%kg/m^3
rhos=2500;           %%kg/m^3
yin=.01;
Tin=600;
k1=6.70*10^10*exp(-12556/T);
K1=65.5*exp(961/T);
r=(0.05*k1*y)/(T*(1+K1*y)^2);
diffy2(1)=(Q*ctot*(yin-y)-alpha*V*r)/(ep*V*ctot);
diffy2(2)=Q*ctot*Mhat*cpg*(Tin-T)+alpha*V*(-delhrxn)*r/(V*(ep*rhog*cpg+(1-ep)*rhos*cps));
diffy2=diffy2';
end

yin=.01;
Tin=600;
zspan=[0 1];
ic=[yin Tin];
[t,y]=ode45(@Ljallo_f_q2,zspan,ic);
plot(t,y(:,1));
xlabel('time');
ylabel('Mol Fraction CO');

Thank you for fixing that for me. I set it from 0 to 1 I'm getting a straight hoorizontal line now. Which is incorrect.

How would I fix this?

[Walter Roberson]

You are not getting a straight horizontal line. You have again made the mistake of not accounting for scale.

Give the command

xlim([0.005 1]); ylim auto

and you will see that the line is not straight at all!