新人求助 用龙库塔迭代求解水面线，程序已编完无法运行

yang

本人最近新学Fortran语言，用龙库塔迭代法求解水面线，程序已编完但无法运行，提示我微分方程设置的语法有问题代码如下：

[Fortran] 纯文本查看 复制代码

查看代码

打印?

01	program main

02	dimension x(0:20),y(0:20),z(0:20),p(0:21),r(0:20),c(0:20)

03	double precision x,y,z,p,h,x0,y0,z0,p0,n,k1,k2,k3,k4,k,q,g,a,b,m,s,r,c

04

k=0.009

05

g=9.8

06

s=0.022

07

m=1.5

08

a=20

09

b=10

10

q=45

11	r=((b+m*y)*y)/(b+2*sqrt(1+m**2)*y)

12	c=(1/s)*r**(1/6)

13	f(x,y)=(k-q**2/(a*c*sqrt(r))**2)/(1-q**2*b/(g*a**3)

14

h=30

15

n=20

16

x(0)=0.0

17

y(0)=1.88

18

r(0)=1.44

19	c(0)=363.6

20	do i=0,n-1

21	r(i+1)=((b+m*y(i+1))*y(i+1))/(b+2*sqrt(1+m**2)*y(i+1))

22	c(i+1)=(1/s)*r(i+1)**(1/6)

23	x(i+1)=(i+1)*h

24	k1=f(x(i),y(i))

25	k2=f(x(i)+0.5*h,y(i)+0.5*h*k1)

26	k3=f(x(i)+0.5*h,y(i)+0.5*h*k2)

27	k4=f(x(i)+h,y(i)+h*k3)

28	y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6

29

end do

30

p(0)=0

31

do i=0,n

32	z(i)=p(i)+exp(-p(i))-1

33	p(i+1)=p(i)+h

34

end do

35	open (10,file="1,txt")

36	write(10,*) (x(i),y(i),z(i),i=0,20)

37

end

li913

13 行， f 未定义。