Mr Singer,

thank you very much for the help, now the script works but above all you have expanded my knowledge on "protection" and "conditions" that it is necessary to impose on functions. I really appreciate it.

Thanks again.

Kind regards,

Marco

Dear Neil,

First of all thanks for the reply and sorry if I cut off the top of the script, my bad. I have now entered the whole script. I probably think that the bug is due to the fact that the rho variable as it is defined does not depend directly on the time. Even seeing your script that has been helpful to me, I am not able to fix it.

g0 = 196133/20000;
R = 287.05;
\[Gamma] = 1.4;
Cpmax = ((\[Gamma] + 1)^2/(4 \[Gamma]))^(\[Gamma]/(\[Gamma] - 
    1))*4/(\[Gamma] + 1);
Rt = 6325766; 

g[z_] := g0*(Rt/(Rt + z))^2;
Cl = 0;
Cd = 1.67929;
S = 12;
m = 5424.9;

z[h_] := r0*h/(r0 + h) /. r0 -> 6325766.
i[z_] := Piecewise[{{1, z <= 11000}, {2, 11000 < z <= 20000}, {3, 
     20000 < z <= 32000}, {4, 32000 < z <= 47000}, {5, 
     47000 < z <= 51000}, {6, 51000 < z <= 80000}, {7, 
     80000 < z < 120000}}];
Tstd[h_] := Module[
  {Lstd, zstd, T0, j, out},

  (*Coefficients definitions*)
  Lstd = {-6.5, 0., 1., 2.8, 0, -2.8, 0};
  zstd = {0., 11000., 20000., 32000., 47000., 51000, 80000};
  T0 = {288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 192.247};

  (*Index i[z]*)
  j = i[z[h]];

  (*ISA model*)
  out = T0[[j]] + Lstd[[j]]/1000 (z[h] - zstd[[j]]);

  Return[out];
  ]
Pstd[h_] := Module[
  {Lstd, zstd, T0, pb, j, out},

  (*Coefficients definitions*)
  Lstd = {-6.5, 0., 1., 2.8, 0, -2.8, 0};
  zstd = {0., 11000., 20000., 32000., 47000., 51000, 80000};
  T0 = {288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 192.247};
  pb = {101325, 22694, 5544.63, 900.192, 120.096, 76.3148, 1.17443};

  (*Index i[z]*)
  j = i[z[h]];

  (*ISA model*)
  out = If[Lstd[[j]] == 0, 
    pb[[j]] Exp[-g0/(R*T0[[j]]) (z[h] - zstd[[j]])], 
    pb[[j]] (1 + Lstd[[j]]/(T0[[j]]*1000) (z[h] - zstd[[j]]))^(-(g0/(
      Lstd[[j]]*0.287)))];
  Return[out];
  ]
\[Rho]std[h_] := Module[
  {out},
  out = Pstd[h]/(R*Tstd[h]);

  Return[out];
  ]

sol = NDSolve[{V'[t] == -0.5* \[Rho]std[h[t]]*(V[t])^2 ((Cd*S)/m) + 
     g[h[t]] Sin[\[Theta][t]],
   \[Theta]'[
     t] == -0.5*\[Rho]std[h[t]]*
      V[t] ((Cl*S)/m) + (g[h[t]]*Cos[\[Theta][t]])/
      V[t] - (V[t]/(Rt + h[t])) Cos[\[Theta][t]],
   h'[t] == -V[t]*Sin[\[Theta][t]],
   V[0] == 11140,
   \[Theta][0] == 6.92*\[Pi]/180,
   h[0] == 120000
   }, {V, \[Theta], h}, {t, 0, 7200}]

Kind regards,

Marco