（#忽略我的用户名。。）我是mathematica新手（从下载到现在才半个月。。）我的作业是模拟弹道轨迹并调初速度和角度来精确打击目标（考虑空气阻力，科氏力，重力的变化），公式都弄出来了，但是代码总是出错。。。
下为代码以及运行结果图片：
ClearAll[y, v, Rx, vx, vy, vz, sx, sy, sz, vx0, vy0, vz0, x, y, z, G, \
alpha, theta, betta, fi, t]
Rx[y_, v_] := 0.5*0.35*v^2*(1 - 2.1904*10^-5*y)^4.4 /; 0 < y <= 9300;
Rx[y_, v_] :=
0.5*0.35*v^2*0.292257*
E^(-2.1206426*(ArcTan[(2.344*(y - 9300) - 6328)/32221.057] +
0.193525)) /; 9300 < y <= 12000;
Rx[y_, v_] :=
0.5*0.35*v^2*0.193725*E^(-(y - 12000)/6483.305)*288.9/221.5 /;
12000 < y <= 30000;
Rx[y_, v_] := 0 /; y > 30000;
G[y_] := 9.8*((6300*10^3)/(6300*10^3 + y))^2 /; y > 0;
Fcorx[vz_, vy_] := -2*
m*\[Omega]*(vz*Sin[fi] + vy*Cos[fi]*Sin[\[Alpha]lpha]);
Fcory[vz_, vx_] := -2*
m*\[Omega]*(-vx*Cos[fi]*Sin[alpha] - vz*Cos[fi]*Cos[alpha]);
Fcorz[vy_, vx_] := -2*m*\[Omega]*(vy*Cos[fi]*Cos[alpha] - vx*Sin[fi]);
m = 3000;
alpha = betta + Pi/2;
\[Omega] = 7.2921*10^-5;
vx0 = v0*Cos[theta]*Sin[betta];
vy0 = v0*Sin[theta];
vz0 = v0*Cos[theta]*Cos[betta];
Manipulate[
Module[{s1 =
NDSolve[{sx''[t] == 1/m*(-Fcorx[vz[t], vy[t]] - Rx[sy[t], vx[t]]),
sy''[t] == 1/m*(-Fcory[vz[t], vx[t]] - Rx[sy[t], vy[t]]),
sz''[t] == 1/m*(-Fcorz[vy[t], vx[t]] - Rx[sy[t], vz[t]]),
sx'[t] == vx[t], sy'[t] == vy[t], sz'[t] == vz[t],
sx'[0] == vx0, vy'[0] == vy0, vz'[0] == vz0, sx[0] == 0,
sy[0] == 0, sz[0] == 0}, {vx, vy, vz, sx, sy, sz}, {t, 0,
300}]}, {ParametricPlot3D[{sx[t], sy[t], sz[t]} /. s1, {t, 0,
"time"},
PlotRange -> {{0, 20000}, {0, 70000}, {0, 70000}}, {{time, 300},
1, 300}]}], {v0, 1200, 3000}, {betta, 0, 2*Pi}, {theta, 0, Pi/2}]

姑且可以这么说：Module只将其内部显式存在的变量局部化，要用就用Block。改完以后还有问题，是方程的问题。
稍微看了一眼，注意到：
Rx[y_, v_] := 0 /; y > 30000;
/;这个东西在解方程的时候少用。要用就用Piecewise（而且最高阶导数不能是Piesewise函数，违背数值积分原理，要做分段就得分段积分，用WhenEvent跳出）