This works (note that Matlab doesn't accept implied multiplication), but doesn't look sensible - what about deriv(2)?
y_initial=[0,101000,0,0]; %[f,smpressure,v,x]
t=0:0.00002:0.02; %sec
% Calling ode45 function
[t,y]=ode45(@(t,y) untitled13(t,y),t,y_initial);
plot(t,y),grid
function [deriv] = untitled13(~,y)
deriv=zeros(4,1);
%Required Data
P=5.000000e-04; % Perforation Diameter (m)
fFT=1.0000; % Force temperature factor
beta=7.5000000e-10; % Burning rate coefficient
alpha=1.0000; % Burning rate exponent
F=1.000000e+06; % Force per unit mass of propellant (J/kg)
C=1.330000e+01; % Mass of propellant (kg)
Vo=2.000000e-02; % Volume of empty chamber (m^3)
A=1.890000e-02; % Area of base of the projectile (m^2)
rho=1600.00; % Density of Propellant (kg/m^3)
b=1.000000e-03; % Covolume of Propellant (m^3/kg)
gamma=1.2600; % Ratio of specific heats for propellant
mp=4.300000e+01; % Mass of projectile (kg)
fR=1.0000; % Downtube resistance factor
Pa=101000; % Pressure in ambient air (0.101 MPa) (Pa)
br=5.000e+06; %Resistive Pressure of the Bore (Pa) [input table]
deriv(1) = -(beta/P)*(y(2)^alpha);
Z = 1-y(1);
Fdot = fFT*F;
y(2) = (C*Z*Fdot-(gamma-1)*(0.5*mp*y(2)))/Vo+(A*y(3)+C*Z*((1/rho)-b)); % C*Z*(... not C*Z(...
Pb = (y(2)+(C*(fR*br))/(3*mp))/(1+(C/(3*mp)));
deriv(3) = (1/mp)*(A*(Pb-Pa)+(fR*br));
deriv(4) = y(3);
end
