I'm trying to solve this code for internal ballistic using ode45.This is my first time using ode

This is the simplified code

Answer 1

[Alan Stevens]

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

[Jayadrath Poojary]

Sorry had to change y(2) to deriv(2). Thank you. It's working now.