моделирование в языковой программе scilab

clear all
m=4;
r(1)=0.015; s=%pi*r(1).^2; g=9.8;
ro=1000; c=0.5; rg=7800;
v0=50;
V=m./rg;
Ft=m*g;
FA=ro*g*V;
k=0.5*c*s*ro;
alpha=50*%pi/180;
vx1(1)=v0*cos(alpha);
vy1(1)=v0*sin(alpha);
ax1(1)=-k*sqrt(vx1(1).^2+vy1(1).^2).*vx1(1)./m;
ay1(1)=g-(k*sqrt(vx1(1).^2+vy1(1).^2).*vy1(1))./m+FA./m;
x1(1)=0;
y1(1)=0;
dt=0.001;
i=1;
for t=0:dt:1.7;
vx1(i+1)=vx1(i)+ax1(i)*dt;
vy1(i+1)=vy1(i)+ay1(i)*dt;
x1(i+1)=x1(i)+vx1(i)*dt+0.5*ax1(i)*dt^2;
y1(i+1)=y1(i)+vy1(i)*dt+0.5*ay1(i)*dt^2;
ax1(i+1)=-k*sqrt(vx1(i).^2+vy1(i).^2).*vx1(i)./m;
ay1(i+1)=g-(k*sqrt(vx1(i).^2+vy1(i).^2).*vy1(i))./m+FA./m;
i=i+1;
end
yx1=max(y1);
xy1=max(x1);
t=0:dt:1.7+dt;
i=1;
x2(1)=0;
y2(1)=0;
dt=0.001
alpha=65*%pi./180;
vx2(1)=v0*cos(alpha);
vy2(1)=v0*sin(alpha);
ax2(1)=-k*sqrt(vx2(1).^2+vy2(1).^2).*vx2(1);
ay2(1)=-(k*sqrt(vx2(1).^2+vy2(1).^2).*vy2(1))./m+FA./m;
for t=0:dt:1.7;
vx2(i+1)=vx2(i)+ax2(i)*dt;
vy2(i+1)=vy2(i)+ay2(i)*dt;
x2(i+1)=x2(i)+vx2(i)*dt+0.5*ax2(i)*dt^2;
y2(i+1)=y2(i)+vy2(i)*dt+0.5*ay2(i)*dt^2;
ax2(i+1)=-k*sqrt(vx2(i).^2+vy2(i).^2).*vx2(i);
ay2(i+1)=g-(k*sqrt(vx2(i).^2+vy2(i).^2).*vy2(i))./m+FA./m;
i=i+1;
end
yx2=max(y2);
xy2=max(x2);
t=0:dt:1.7+dt;
i=1;
x3(1)=0;
y3(1)=0;
dt=0.001
alpha=70*%pi/180;
vx3(1)=v0*cos(alpha);
vy3(1)=v0*sin(alpha);
ax3(1)=-k*sqrt(vx3(1).^2+vy3(1).^2).*vx3(1)./m;
ay3(1)=g-(k*sqrt(vx3(1).^2+vy3(1).^2).*vy3(1))./m+FA./m;
for t=0:dt:1.7;
vx3(i+1)=vx3(i)+ax3(i)*dt;
vy3(i+1)=vy3(i)+ay3(i)*dt;
x3(i+1)=x3(i)+vx3(i)*dt+0.5*ax3(i)*dt^2;
y3(i+1)=y3(i)+vy3(i)*dt+0.5*ay3(i)*dt^2;
ax3(i+1)=-k*sqrt(vx3(i).^2+vy3(i).^2).*vx3(i)./m;
ay3(i+1)=g-(k*sqrt(vx3(i).^2+vy3(i).^2).*vy3(i))./m+FA./m;
i=i+1;
end
yx3=max(y3);
xy3=max(x3);
t=0:dt:1.7+dt;
xy=[xy1,xy2,xy3];
yx=[yx1,yx2,yx3];
figure(1)
subplot(2,2,1)
plot(t,ax1,t,ax2,t,ax3);
xgrid
xlabel('t');
ylabel('ax1');

subplot(2,2,2)
plot(t,ay1,t,ay2,t,ay3);
xgrid
xlabel('t');
ylabel('ay1');

subplot(2,2,3)
plot(t,x1,t,x2,t,x3);
xgrid
xlabel('t');
ylabel('x');

subplot(2,2,4)
plot(t,y1,t,y2,t,y3);
xgrid
xlabel('t');
ylabel('y');

figure(2)
subplot(2,2,1)
plot(t,vx1,t,vx2,t,vx3);
xgrid
xlabel('t');
ylabel('vx');

subplot(2,2,2)
plot(t,vy1,t,vy2,t,vy3);
xgrid
xlabel('t');
ylabel('vy');
есть недоработанный код.нужно сделать 20 вариант