the plot of start is good but i do not know why the other plot doesn't want to plot steady state +final value. I hope someone will fix it thanks in advance

R=8.314/32;

T0=2930;

a=(12)/((2.027*10^6)^0.45);

rhoP=1920;

Astar=pi*0.25^2;

k=1.35;

n=0.45;

P0=101325;

P = P0;

%syms P(t)

%at beginning of the integration set initial values for the persistent variables

rp=0.35; %initial port radius

t1=0; %initial time step

rhoO=P/(R*T0); %gas density

rp=min(rp+((a*P^n)*10^-3)*(t1-t1),0.7);

Ab=2*pi*rp*8;%burn area

v0 = 0.1;

dP=@(t,P)Fun(t,P,R,T0,rp,a,n,t1,Ab,rhoP,Astar,k,v0);%@(t,P)(Ab.*a.*P.^n.*(rhoP-rhoO)-P.*Astar.*sqrt(k/(R.*T0)).*(2/(k+1)).^((k+1)/(2.*(k-1)))).*R.*T0./v0;

[t,P]=ode45(dP, [0,0.1], P);

y = P;

figure(1)

plot(t,y)

xlabel("Time (s)")

ylabel("Chamber Pressure (Pa)")

title("Chamber Pressure vs Time (Start-Up)")

dP=@(t,P)Fun(t,P,R,T0,rp,a,n,t1,Ab,rhoP,Astar,k,v0);%@(t,P)(Ab.*a.*P.^n.*(rhoP-rhoO)-P.*Astar.*sqrt(k/(R.*T0)).*(2/(k+1)).^((k+1)/(2.*(k-1)))).*R.*T0./v0;

[t,P]=ode45(dP, [0.1,60], P);

y = P;

figure(2)

plot(t,y)

xlabel("Time (s)")

ylabel("Chamber Pressure (Pa)")

title("Chamber Pressure vs Time ")

function dP = Fun(t,P,R,T0,rp,a,n,t1,Ab,rhoP,Astar,k,v0)

char t1;char rp;

if t==0

rp=0.35;

t1=0;

end

if t>=0.1

rp>=0.7;

Ab=0;

end

v0=pi*rp^2*8; %recalculating free volume available in chamber

t1=t; %storing current time step to be used in next iteration

rhoO=P/(R*T0); %gas density

rp=min(rp+((a*P^n)*10^-3)*(t-t1),0.7);

Ab=2*pi*rp*8;%burn area

dP = (Ab.*a.*P.^n.*(rhoP-rhoO)-P.*Astar.*sqrt(k/(R.*T0)).*(2/(k+1)).^((k+1)/(2.*(k-1)))).*R.*T0./v0;

end

[Sai Sri Pathuri]

When the script is calling dp function for second time, P is a matrix. Hence the equation

rp=min(rp+((a*P^n)*10^-3)*(t-t1),0.7);  

is to be replaced with below equation to support calculation of power of a matrix (P.^n)

rp=min(rp+((a*P.^n)*10^-3)*(t-t1),0.7);

Then second graph can be plotted