Solving nonlinear second order diferential equation

Dados iniciais

clear
Massa_foguete_seco        = 3000; % Kg
Massa_combustivel_inicial = 3000; % kg
Empuxo                    = 300000; % Newton (sentido de y)
Consumo_combustivel       = 100; % Kg/segundo
C_arrasto                 = 50; % Newton/(metros/segundo)
Gravidade                 = 9.81; % Metros/segundo^2 (sentido contrario de y)
rho                       = 2; % Kg/metros cubicos
Area_secao                = 15; % Metros quadrados
                              

Peso do Foguete.

syms y(t) t

Massa_combustivel = Massa_combustivel_inicial - Consumo_combustivel * t;

Massa_foguete = Massa_combustivel + Massa_foguete_seco;

Peso_foguete = Massa_foguete * -Gravidade;

Cálculo da Aceleração

Forca_arrasto    = -0.5 * C_arrasto * rho * Area_secao * diff(y,t)^2;
Forca_empuxo     = Empuxo;
Peso             = Peso_foguete;

Forca_resultante = Forca_arrasto + Forca_empuxo + Peso;

Equações do movimento

Aceleracao = Forca_resultante/ Massa_foguete

ode = diff(y,t,2) == Aceleracao;

dy = diff(y);
cond1 = y(0)  == 0;
cond2 = dy(0) == 0;
cond = [cond1 cond2];

Altitude    = dsolve(ode,cond)
Velocidade = diff(Altitude);
Aceleracao = diff(Altitude,2);

Answer 1

[Stephan]

%% Dados iniciais
clear
Massa_foguete_seco        = 3000; % Kg
Massa_combustivel_inicial = 3000; % kg
Empuxo                    = 300000; % Newton (sentido de y)
Consumo_combustivel       = 100; % Kg/segundo
C_arrasto                 = 50; % Newton/(metros/segundo)
Gravidade                 = 9.81; % Metros/segundo^2 (sentido contrario de y)
rho                       = 2; % Kg/metros cubicos
Area_secao                = 15; % Metros quadrados

%% Peso do Foguete.
syms y(t)
Massa_combustivel = Massa_combustivel_inicial - Consumo_combustivel * t;
Massa_foguete = Massa_combustivel + Massa_foguete_seco;
Peso_foguete = Massa_foguete * -Gravidade;

%% Cálculo da Aceleração
Forca_arrasto    = -0.5 * C_arrasto * rho * Area_secao * diff(y,t)^2;
Forca_empuxo     = Empuxo;
Peso             = Peso_foguete;
Forca_resultante = Forca_arrasto + Forca_empuxo + Peso;

%% Equações do movimento
Aceleracao = Forca_resultante/ Massa_foguete;
ode = diff(y,t,2) == Aceleracao

%% Solve numerically
[eqs, vars] = odeToVectorField(ode);
odefun = matlabFunction(eqs,'vars',{'t','Y'})

tspan = [0 5];
init = [0 0];
[t,y] = ode45(odefun,tspan,init);
 
plot(t,y)

[Rafael Pessoa]

THANK U, it worked, but how can i derivate this function? and why there are two lines in the graphic?

[Stephan]

This is a numerical solution. One line represents the velocity, the other the way. Since this is a second order system, which is transformed into a first order system of 2 ode, you get the solution to both, the way and the velocity.

See here also:

https://www.wolframalpha.com/input/?i=y%27%27%28t%29+%3D+-%28981*t+-+750*y%27%28t%29%5E2+%2B+241140%29%2F%28100*t+-+6000%29%2C+y%280%29+%3D+0%2C+y%27%280%29%3D0

To calculate the acceleration from the results use:

%% Dados iniciais
clear
Massa_foguete_seco        = 3000; % Kg
Massa_combustivel_inicial = 3000; % kg
Empuxo                    = 300000; % Newton (sentido de y)
Consumo_combustivel       = 100; % Kg/segundo
C_arrasto                 = 50; % Newton/(metros/segundo)
Gravidade                 = 9.81; % Metros/segundo^2 (sentido contrario de y)
rho                       = 2; % Kg/metros cubicos
Area_secao                = 15; % Metros quadrados

%% Peso do Foguete.
syms y(t)
Massa_combustivel = Massa_combustivel_inicial - Consumo_combustivel * t;
Massa_foguete = Massa_combustivel + Massa_foguete_seco;
Peso_foguete = Massa_foguete * -Gravidade;

%% Cálculo da Aceleração
Forca_arrasto    = -0.5 * C_arrasto * rho * Area_secao * diff(y,t)^2;
Forca_empuxo     = Empuxo;
Peso             = Peso_foguete;
Forca_resultante = Forca_arrasto + Forca_empuxo + Peso;

%% Equações do movimento
Aceleracao_num = Forca_resultante/ Massa_foguete;
ode = diff(y,t,2) == Aceleracao_num;

%% Solve numerically
[eqs, vars] = odeToVectorField(ode);
odefun = matlabFunction(eqs,'vars',{'t','Y'});

tspan = [0 5];
init = [0 0];
[t,y] = ode45(odefun,tspan,init);

% Calculate acceleration from numeric solution - Create ab function handle
% that represents the acceleration ode, but now we know the numeric
% solution for the ode system, so we can use the velocity informations
% needed in the original ode.
Aceleracao_num = @(t,Y2) -(981*t - 750*Y2.^2 + 241140)./(100*t - 6000);

% Calculate those values based on the results from ode45 and also save them
% in y-vector, column 3
y(:,3) = Aceleracao_num(t,y(:,2));

% plot the whole system, with y, y_dot and y_doubledot
plot(t,y)