[求助]我编写的导弹拦截微分方程组出现了问题

[深蓝浅韵]

matlab的各位牛人，走过路过的看一下我的问题
我编写的m文件建立弹道导弹拦截模型出现了很多问题，有没有谁看出来，帮我修改一下
我是matlab的初学者，马上要交论文了，现在是火烧眉毛啊
多谢各位了！！！
首先是微分方程的函数定义
function drdqdzdadt=ydt(t,Y);
global at A b V Vt K r q z a;
Y=[r,q,z,a];
A=atan(b*z/r);
drdqdzdadt=[Vt*cos(q-at)-V*cos(A)*cos(q-a);(V*cos(A)*sin(q-a)-Vt*sin(q-at))/r;-V*sin(A);K*(V*cos(A)*sin(q-a)-Vt*sin(q-at))/r]
%{
Vt：目标速率；
V：导弹速率；
K：攻击平面导引比例；
b：攻击平面法向导引比例；
A：θ；
a:σ;
%}
然后是主程序
Vt=input('请输入目标速率');
at=input('请输入目标速度与水平面的夹角');
V=1000;
b=10;
K=1;
x(1)=500;
y(1)=500;
[t,Y]=ode45(@ydt,[0:0.1:30],[15000;pi/4;3000;pi/6]);
ptr=[x(1)-Vt*cos(at)*t;y(1)-Vt*sin(at)*t;z(1)];
pmr=[ptr(1,-r*cos(q);ptr(2,-r*sin(q);z(1)-r*tan(A)];
D=det(ptr-pmr);
output=[t' r' z' D' ptr(2,']
建模的方程组在附件中
r——目标线（导弹与目标的连线）和基线（目标线在oxz平面内的投影）

[anbcjys]

cos(q-at)是什么 另外没这么解过方程 我是用simulink做的导弹拦截