I am doing a project on the motion of spherical objects through air. It has 3 equations which need to be solved using the ODE solver. The SS of Matlab Code is provided. In the Matlab Code equations have variables that are matrices. When i try to replicate that in Scilab it gives me error.

Matrix need to be used as equations have the same variable "velocity" with different values in x,y and z directions. All depending on the variable against which we have taken derivative and trying to solve the equation. This is distance the independent variable

Please Note that the rest of the code is FIne in Scilab as it works when i change the matrix to a single variable. However just tp be sure i have provided the complete scilab code below:

Can Anyone Please Help me

Sci Lab Code:

clear all
clc

D=input('What is the diameter of the ball under consideration ')

vx=input("What is the initail velocity in x-axis ")
vy=input("What is the initail velocity in y-axis ")
vz=input("What is the initail velocity in z-axis ")

vvector=[vx, vy, vz];
vmag=sqrt(vx^2+vy^2+vz^2);

sx=input('What is the initial spin rate in x-axis ')
sy=input('What is the initial spin rate in y-axis ')
sz=input('What is the initial spin rate in z-axis ')

smag=sqrt(sx^2+sy^2+sz^2);
svector=[sx, sy, sz];

C_d=0.45;
C_l=3.19*10^-1*(1-exp(-2.48*10^3*smag));
g=9.81;

wx=input('What is the value of wind in x-axis ')
wy=input('What is the value of wind in y-axis ')
wz=input('What is the value of wind in z-axis ')

wvector=[wx, wy, wz];
wmag=sqrt(wx^2+wy^2+wz^2);

// These are the constants that will be used in the equations

roh=1.22;
Area=(%pi*D^2)/4;

vwresultant=vvector-wvector;
vwmag=sqrt(vwresultant(1)^2+vwresultant(2)^2+vwresultant(3)^2);

//Here we will define the equations that will solve for the velocity in x, y and z direction

dy=zeros(6,1);
v=vvector;

funcprot(0)
function dy(2)=fvx(t1,v)
   dy(1)=v(1);
   dy(2)=((-1/2)*roh*Area*vwmag)*(C_d*(v(1)-wx)-C_l*(sy*(v(z)-wz)-sz*(v(2)-wy))/smag);
endfunction

function dy(4)=fvy(t2,v)
    dy(3)=v(2);
    dy(4)=((-1/2)*roh*Area*vwmag)*(C_d*(v(2)-wy)-C_l*(sz*(v(1)-wx)-sx*(v(3)-wz))/smag);
endfunction

function dy(6)=fvz(t3,v)
    dy(5)=v(3);
    Velocityz=((-1/2)*roh*Area*vwmag)*(C_d*(v(3)-wz)-C_l*(sx*(v(2)-wy)-sy*(v(1)-wx))/smag)-g;
endfunction

//Now we will define the initial conditions and then solve
//Uptil this point we have simply prompt for values of known conditions, and plugged in equations of motion as functions. Now what we need to do is to integrate to find the values of velocity at particular points and then use them to find the velocity at next points.

vvaluesx(1)=vx;
vvaluesy(1)=vy;
vvaluesz(1)=vz;

for rangev=1:1:100 
    vmagvalues(rangev)=sqrt(vvaluesx(rangev)^2+vvaluesy(rangev)^2+vvaluesz(rangev)^2);
    
    vvaluesx(rangev+1)=ode(vvaluesx(rangev),((rangev-1)/10),(rangev/10),fvx);
    vvaluesy(rangev+1)=ode(vvaluesy(rangev),((rangev-1)/10),(rangev/10),fvy);
    vvaluesz(rangev+1)=ode(vvaluesz(rangev),((rangev-1)/10),(rangev/10),fvz);
    
    vx=vvaluesx(rangev+1);
    vy=vvaluesy(rangev+1);
    vz=vvaluesz(rangev+1);
end

//Now we will plot the results to see the trend

t=0:0.1:9.9; //We have started from 0 to get equal values since above we get 100+1 values of vvalues variables. We want the values to be equal
plot(vvaluesx);
plot(vvaluesy);
plot(vvaluesz);