对于Ancona的鲨鱼与鱼的微分方程，求任意T时刻的x1,x2
方程已上传   下面贴上我写出来的程序

subroutine gelr(t,y,m,h,n,z,f,d)         !这里是用改进的欧拉公式，用定步长进行全区间积分
dimension y(m),z(m,n),d(m)               !t是微分方程组的初值点，即积分的起点
double precision y,z,d,t,h,x             !y是一维数组，长度为M，微分方程组在初始点T出的M歌未知函数的
do i=1,m                                 !初值，返回值在z(i,1)中
z(i,1)=y(i)                              !H，从初值点开始，对微分方程积分一步的步长
end do                                   !N，对微分方程积分的步数
do 50 j=2,n                              !Z，二维数组，体积为M*N，存放N个积分点上的未知函数值  
   x=t+(j-2)*h                           !F，子程序，用于计算方程组中各方程右端函数值
   call f(x,y,m,d)                       !D，一维数组，长度为M，存放右端函数值
   do 20 i=1,m
  20    y(i)=z(i,j-1)+h*d(i)
      x=t+(j-1)*h
      call f(x,y,m,d)
          do 30 i=1,m
         30    d(i)=z(i,j-1)+h*d(i)
                do 40 i=1,m
                y(i)=(y(i)+d(i))/2.0
                z(i,j)=y(i)
                
           
    40 continue
50 continue
return
end
external f
dimension y(2),z(2,50),d(2)
double precision y,z,d,t,h,x             !这里是我计算是用的一些初值
t=0
y(1)=25.0
y(2)=2.0 
h=0.1
m=2
n=50
call gelr(t,y,m,h,n,z,f,d)
do i=1,n
   write(*,30)(z(j,i),j=1,m)
end do
30 format(1x,'y(1)=',f14.1,3x,'y(2)=',f4.1)
pause
end

subroutine f(t,y,m,d)              !这里是对于dx1/dt,dx2/dt的定义，c1,c2对应方程里面的λ1,λ2
dimension y(2),d(2),r(2),c(2)
double precision y,d
r(1)=1
r(2)=0.5
c(1)=0.1
c(2)=0.02
d(1)=y(1)*(r(1)-c(1)*y(2))
d(2)=y(2)*(-r(2)+c(2)*y(1))
return
end



这个程序运行后。其实计算结果是对的，但是要求N不能是50，必须要更大
问题就出在这里，大概50之后，这个程序计算的结果就是x1<0,x2很小了
这个与正确答案是有很大偏差的，
所以想请有经验的大大，帮我看看，是哪个方面除了问题
我是初学者，自己很难查出有什么错误来