function f=fun1(time,y)
vx=8.0;vy=8.0;vz=0;
ax=0.0;ay=0.0;az=0.0;
latitude=45*pi/180;
R=6371020.0;
wie=7.2921158*10^(-5);
g=9.80;
fx=ax-(2*wie*sin(latitude)+vx/R*tan(latitude))*vy;
fy=ay+(2*wie*sin(latitude)+vx/R*tan(latitude))*vx;
fz=az-vy^2/R-(2*wie*cos(latitude)+vx/R)*vx+g;
a12=wie*sin(latitude)+vx/R*tan(latitude);
a13=-(wie*cos(latitude)+vx/R);
a15=-1/R;
a21=-a12;
a23=-vy/R;
a24=-a15;
a26=-wie*sin(latitude);
a31=-a13;
a32=-a23;
a34=tan(latitude)/R;
a36=wie*cos(latitude);
a42=-fz;
a43=fy;
a44=vy/R*tan(latitude)-vz/R;
a45=2*wie*sin(latitude)+vx/R*tan(latitude);
a46=(2*wie*cos(latitude)+vx/R*(sec(latitude))^2)*vy+2*wie*sin(latitude)*vz;
a51=-a42;
a53=-fx;
a54=-2*(wie*sin(latitude)+tan(latitude)/R);
a55=-vz/R;
a56=-(2*wie*cos(latitude)+vx/R*(sec(latitude))^2)*vx;
a65=1/R;
A=[0 a12 a13 0 a15 0;
   a21 0 a23 a24 0 a26;
   a31 a32 0 a34 0 a36;
   0 a42 a43 a44 a45 a46;
   a51 0 a53 a54 a55 a56;
   0 0 0 0 a65 0];
C=[-0.01*pi/180/3600;-0.01*pi/180/3600;-0.01*pi/180/3600;10^(-4)*g;10^(-4)*g;0];
f=A*y+C;