module functions
implicit none
real,parameter:: beta = 0.001, ni = 0.1
contains
function func(x,y)
real, intent(in)::x
real,dimension(:),intent(in):: y
real,dimension(size(y))::func
func(1)= -beta*y(2)*y(1)
func(2)= beta*y(1)*y(2) - ni*y(2)
func(3)= ni*y(2)
end function func
end  module functions

program sir
use functions
implicit none
real :: h,x,diff
real, dimension(3) :: y,k1,k2,k3,k4,ynew
integer :: i,n
write (*,*) 'h= '
read (*,*) h
write (*,*) 'n= '
read (*,*) n
write (*,*) 'S init = '
read(*,*) y(1)
write(*,*) 'I init = '
read(*,*) y(2)
write(*,*) 'R init = '
read (*,*) y(3)
x=0.e0
write(*,*) x,y(1),y(2),y(3)
do i=0,n-1
x = x+h
k1(:)=h*func(x,y)
k2(:)=h*func(x+h/2,y+k1/2)
k3(:)=h*func(x+h/2,y+ k2/2)
k4(:)=h*func(x+h,y+k3)
ynew(:)= y (:)+ (k1 +2*k2+2* k3+k4)/6
y(:)= ynew (:)
write(*,*) x,y(1),y(2),y(3)
end do
end program