program main

    parameter(pi=3.1415926,e=2.1e9,r=0.01,dx=0.1,p=7.85e9,dt=0.15)
    real k(10,10),m(10,10),b(10,10),f(25,10),u(27,10)
    real a1,b1
    integer a(2,2)
    integer i,j
    real area,gxs,zxs,g1,g2

c波速、总长度、横截面积
    c=sqrt(e/p)
    l=(10-1)*dx/c+2
    area=pi*(r/2)*(r/2)

    gxs=(e*area)/dx
    zxs=(p*area*dx)/2
    g1=gxs
    g2=(-1)*gxs

      data((k(i,j), i=1,10),j=1,10) /100*0/
      a(1,1)=g1
    a(1,2)=g2
    a(2,1)=g2
    a(2,2)=g1
      do i=1,10
     k(i,i)=k(i,i)+a(1,1)
      
     enddo
     do j=2,9
        k(j,j)=k(j,j)+a(2,2)
       enddo
    do i=1,9
      k(i,i+1)=a(1,2)
      k(i+1,i)=a(1,2)
     enddo

    data((m(i,j), i=1,10),j=1,10) /100*0/
    do i=1,10
     m(i,i)=m(i,i)+zxs
      
     enddo


      a1=2/zxs
      b1=(dt*dt)/zxs      
        b=a1*m-b1*k
    data((f(j,i), i=2,10),j=1,25) /225*0/


            do j=1,25
               do i=1,10
       u(0,i)=0
       u(1,i)=0
    
      f(j,1)=y(j*dt)

    u(j+1,i)=b(i,i)*u(j,i)+b(i,i+1)*u(j,i+1)-u(j-1,i)+f(j,i)

    print*,u(j+1,i)   
           enddo
       enddo
    
    

    end
中间的子函数y(x)的程序为：
    function y(x)
    real x,pi,s
    s=1.0e4
    pi=3.1415926
    if((x.gt.0).and.(x.le.2))then
    y=s*sin(pi*x/2)
    else
    y=0
    endif
    end
求各位帮忙！！！