Trying to use an ODE solver for a convergent divergent nozzle. I don't understand why it's telling me the number of elements on the left and right side are different and its not showing anything in the workspace for me to figure it out.

Sorry I'm not one bit talented in coding and there may be a lot of unnecessary lines of code but here goes..

This is the function to be called:

function dxdt=tobecalled(t,x)
pi=3.142;
%Area change:
if 0<t<5
    A=1.75-(0.75*cos(((0.2*t)-1)*pi));    
else
    A=1.25-(0.25*cos(((0.2*t)-1)*pi));
end
D=sqrt(A*(4/pi));    %Diameter formula
g=1.4;  %Gamma
f=0.005;    %Wall friction parameter
T0=293;     %Temperature (adiabatic and constant)
mdot=0.334; %Mass flow rate
mSquare=((x(1))^2);    %Mach square 
phi=1+(((g-1)/2)*mSquare); %Phi constant
cM=(phi*mSquare)/(1-mSquare);   %dM/dx constant
cP=(g*mSquare*x(2))/(1-mSquare);    %dP/dx constant
cRho=(mSquare*x(3))/(1-mSquare);    %dRho/dx constant
cT=((g-1)*mSquare*x(4))/(1-mSquare);    %dT/dx constant
cV=(x(5))/(1-mSquare);  %dV/dx constant
dA=diff(A);   %dAdx
dmdot=diff(mdot); %dmdotdx
dT0=diff(T0);   %dT0/dx
fric=(2*f)/D;   %Friction term
dxdt=zeros(5,1);
dxdt(1)=cM*(((-1/A)*dA)+((g*mSquare)*fric)+(((1+g)/(2*T0))*dT0)+(((1+(g*mSquare))/mdot)*dmdot));    %dM/dx
dxdt(2)=cP*(((1/A)*dA)-((1+((g-1)*mSquare))*(fric))-((phi/T0)*dT0)-(((2*phi)/mdot)*dmdot)); %dP/dx
dxdt(3)=cRho*(((1/A)*dA)-(g*fric)-((phi/(mSquare*T0))*dT0)-(((g+1)/mdot)*dmdot));    %dRho/dx
dxdt(4)=cT*(((1/A)*dA)-((g*mSquare)*fric)-(((1+(g*mSquare))/mdot)*dmdot))+((((1-(g*mSquare))*phi*x(4))/(1-(mSquare*T0)))*dT0);     %dT/dx
dxdt(5)=cV*(((1/A)*dA)+((g*mSquare)*fric)+((phi/T0)*dT0)+(((1+(g*mSquare))/mdot)*dmdot));   %dV/dx
return

And this is the function which calls the one above and to plot:

close all
[t,x]=ode45(@tobecalled, [0,10],[0 0 0 0 0]');
pi=3.142;
if 0<t<5
    A=1.75-(0.75*cos(((0.2*t)-1)*pi));
else
    A=1.25-(0.25*cos(((0.2*t)-1)*pi));
end
D=sqrt(A*(4/pi));    %Diameter formula
g=1.4;  %Gamma
f=0.005;    %Wall friction parameter
T0=293;     %Temperature (adiabatic and constant)
mdot=0.334; %Mass flow rate
mSquare=((x(1))^2);    %Mach square constant
phi=1+(((g-1)/2)*mSquare); %Phi constant
cM=(phi*mSquare)/(1-mSquare);   %dM/dx constant
cP=(g*mSquare*x(2))/(1-mSquare);    %dP/dx constant
cRho=(mSquare*x(3))/(1-mSquare);    %dRho/dx constant
cT=((g-1)*mSquare*x(4))/(1-mSquare);    %dT/dx constant
cV=(x(5))/(1-mSquare);  %dV/dx constant
dA=diff(A);   %dAdx
dmdot=diff(mdot); %dmdotdx
dT0=diff(T0);   %dT0/dx
fric=(2*f)/D;   %Friction term
plot(t,x(:,1),'LineWidth',0.5); grid on ; hold on;

Answer 1

[darova]

Create breakpoint before the line you have an error using F12. You will find out that variables you are trying to use are empty

This happens because of these lines

dA=diff(A);   %dAdx
dmdot=diff(mdot); %dmdotdx
dT0=diff(T0);   %dT0/dx

You are trying to calculate difference of constant value

T0=293;     %Temperature (adiabatic and constant)
mdot=0.334; %Mass flow rate

Here is an example

diff(2)
ans =
     []