SUBROUTINE difeq(k,k1,k2,jsf,is1,isf,indexv,ne,s,nsi,nsj,y,nyj, *nyk) INTEGER is1,isf,jsf,k,k1,k2,ne,nsi,nsj,nyj,nyk,indexv(nyj),M REAL s(nsi,nsj),y(nyj,nyk) COMMON /sfrcom/ x,h,mm,n,c2,anorm PARAMETER (M=41) INTEGER mm,n REAL anorm,c2,h,temp,temp2,x(M) if(k.eq.k1) then if(mod(n+mm,2).eq.1)then s(3,3+indexv(1))=1. s(3,3+indexv(2))=0. s(3,3+indexv(3))=0. s(3,jsf)=y(1,1) else s(3,3+indexv(1))=0. s(3,3+indexv(2))=1. s(3,3+indexv(3))=0. s(3,jsf)=y(2,1) endif else if(k.gt.k2) then s(1,3+indexv(1))=-(y(3,M)-c2)/(2.*(mm+1.)) s(1,3+indexv(2))=1. s(1,3+indexv(3))=-y(1,M)/(2.*(mm+1.)) s(1,jsf)=y(2,M)-(y(3,M)-c2)*y(1,M)/(2.*(mm+1.)) s(2,3+indexv(1))=1. s(2,3+indexv(2))=0. s(2,3+indexv(3))=0. s(2,jsf)=y(1,M)-anorm else s(1,indexv(1))=-1. s(1,indexv(2))=-.5*h s(1,indexv(3))=0. s(1,3+indexv(1))=1. s(1,3+indexv(2))=-.5*h s(1,3+indexv(3))=0. temp=h/(1.-(x(k)+x(k-1))**2*.25) temp2=.5*(y(3,k)+y(3,k-1))-c2*.25*(x(k)+x(k-1))**2 s(2,indexv(1))=temp*temp2*.5 s(2,indexv(2))=-1.-.5*temp*(mm+1.)*(x(k)+x(k-1)) s(2,indexv(3))=.25*temp*(y(1,k)+y(1,k-1)) s(2,3+indexv(1))=s(2,indexv(1)) s(2,3+indexv(2))=2.+s(2,indexv(2)) s(2,3+indexv(3))=s(2,indexv(3)) s(3,indexv(1))=0. s(3,indexv(2))=0. s(3,indexv(3))=-1. s(3,3+indexv(1))=0. s(3,3+indexv(2))=0. s(3,3+indexv(3))=1. s(1,jsf)=y(1,k)-y(1,k-1)-.5*h*(y(2,k)+y(2,k-1)) s(2,jsf)=y(2,k)-y(2,k-1)-temp*((x(k)+x(k-1))*.5*(mm+1.)*(y(2,k)+ *y(2,k-1))-temp2*.5*(y(1,k)+y(1,k-1))) s(3,jsf)=y(3,k)-y(3,k-1) endif return END