ode45 doesn't run?

I am doing a system of equations that are being solved by ode45. I have six state variables but whenever I try to run the ode45, I get an answer where my x array only contains the initial conditions!! I need to know why doesn't ode45 do over the equations by adding the regular time steps.

Here is the whole code and I apologize for putting it in its entirety:

%Pneumatic Launcher Molecules

if true
  % code
end
clear; clc; close all
%Constants 
global Gg T K Z Cv B Vo D L_tot m A gama_max l h d Patm Ptank Q Pt Pb 
Gg = 1;         %Specific gravity of air
T  = 293;       %Temperature, k
K  = 1.380e-23; %Boltzmann Constant, J/K
Z  = 1;         %Compressiblity Factor
Cv = 1.93;      %Flow Coefficient
B  = 3.11e19;   %Engineering Flow Constant converts quantity into units of molecules per time, sqrt(k)/(Pa*s)
Vo = 0.004196;  %Tank volume, m^3
D  = 0.01913;   %Barrel Diameter, m
m  = 0.0194;    %Mass of the payload, kg
A  = pi*D^2/4;  %Barrel cross sectional area, m^2
gama_max = 0.8; %Flow Ratio (For chocked vs unchocked)
L_tot = 0.916;  %total length, m 
h = 0.04802;    %Slug Length, m
l = 0.8825;     %initial length through the barrel, m  
d = L_tot - l;  %Initial distance from the valve to the slug, m 

Patm = 101300;  %Initial barrel pressure, atmospheric, Pa
Ptank = 300000;             %Ininitial Tank Pressure, Pa
%Initial Conditions

x_ini    = 0;                  %initial distance
xdot_ini = 0;              %initial velocity
Pb_ini   = Patm;            %initial barrel pressure
Pt_ini   = Ptank;           %initial tank pressure
N_ini    = Ptank*Vo/(K*T);  %initial number of molecules in the tank
Nb_ini   = Patm*A*d/(K*T);  %initial number of molecules in the barrel

ini      = [x_ini;  xdot_ini; N_ini; Nb_ini; Pt_ini; Pb_ini]; %Initial Conditions vector

%The modeling equation

Tspan = 0:0.01:0.5;  
options=odeset('event',@stopsim);
[t,x] = ode45(@Plauncher, Tspan, ini, options);

function [dx] = Plauncher(t,x)

global Gg T K Z Cv B Vo D L_tot m A gama_max l h d Patm Ptank Q Pt Pb 

%This function takes the state variables and outputs the required
%parameters
ssure, Pa
%STATE VARIABLES

dx  =zeros(6,1);
X   =x(1); %Position
xdot=x(2); %Velocity
N   =x(3); %Number of Molecules in the tank
Nb  =x(4); %Number of Molecules in the barrel
Pt  =x(5); %Tank Pressure
Pb  =x(6); %Barrel Pressure

gama1 = (Pt-Pb)/Pt;    
if gama1 >= gama_max     
    gama1 = gama_max;    
else     
    gama1 = (Pt-Pb)/Pt;     
end

if Pt < Pb/(1-gama_max)   
    Q = B*Pt*Cv*(1-gama1/(3*gama_max))*sqrt(gama1/(Gg*K*T));   
else   
    Q = (2/3)*B*Pt*Cv*sqrt(gama_max/(Gg*K*T));   
end

dx(1) = xdot;  % Velocity
dx(2) = (A*Pb - A*Patm)/m;  % Acceleration
dx(3) = -Q; %Tank Molecules number Differential
dx(4) =  Q; %Barrel Molecules Number Differential 
dx(5) = -Q*K*T/Vo; % Tank Pressure Differential  
dx(6) =  Q*K*T/(A*(d+X+xdot)); %Barrel Pressure Differential

end

function [ Val, Ister, Dir ] = stopsim( t,x )
%STOPSIM Summary of this function goes here
%   Detailed explanation goes here

global l Pt Pb 

Val(1)=x(1)-l; %Stops whenever the slug is at the exit of the barrel
Val(1)=Pt-Pb; %Stops whenever the tank pressure is less than barrel pressure
Ister(1)=1;
Dir(1)=1;

end

finally, here is what I get at the end:

Warning: Failure at t=3.789533e-13. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (8.077936e-28) at time t. > In ode45 (line 308) In pneumatics_molecules (line 42)

Answer 1

[Prashant Arora]

I understand that you are trying to solve a non-linear ode with some stop intergration constraints defined in ‘stopsim’ function as an ‘event’ property.

To answer your first question, the ordinary differential equation solver functions provided with MATLAB employ a variety of variable-step methods. To get a fixed step solver, you can check the following link.:

https://www.mathworks.com/matlabcentral/answers/98293-is-there-a-fixed-step-ordinary-differential-equation-ode-solver-in-matlab-8-0-r2012b

Regarding the error you mentioned at the end:

I executed your code and realized that ‘Pb’ , ‘Pt’ and ‘Q’ were obtaining imaginary values at the end of execution. It meant that ‘gama1’ obtained negative values, which were later passed on to the “sqrt” function becoming a source of imaginary numbers. This should not have happened as you tried to control this phenomenon using the ‘stopsim’ function. In the 'stopsim' function, you defined the ‘Dir’ to be 1, which means you are tracking only zeros resulting due to an increasing value of event. However, the value of ‘Pt-Pb’ is a decreasing and does not get tracked. Also, the ‘Val’ variable is overwritten in the function. I suggest changing the ‘stopsim’ function to:

if round(x(1)-l)==0 || Pt<Pb
Val(1) = 0;
else
Val(1) = 1;
end
Ister(1) = 1;
Dir(1) = 0; %Track all zeros

The simulation should run now, but will terminate according to the conditions defined in ‘stopsim’, which happens to be at a very early stage of simulation. I would recommend checking the initial conditions and probably equations and units also, to prevent the ‘stopsim’ conditions resulting true very soon.