c DE with variable coefficients
c Differential Equations Ref. Morris & Brown page 94
c x**2(d2y/dx2)-3.*x(dy/dx) -5.*y=0
data c1,c2,xi,xf,nstep/1.,1.,1.,2.,2000/
c one solution is ysol(x)= c1*x**5 +c2/x
ysol(x)= c1*x**5 +c2/x
dx=(xf-xi)/float(nstep)
c initial values are taken at x=dx
y0=ysol(xi)
dydx=5.*c1*xi**4 -c2/xi**2
y1=y0+dx*dydx
kp=int(float(nstep)/40.)
kount=kp
print*,'nstep,dx=',nstep,dx
print*,' '
Print*,' x ynum ysol '
print 100,xi,y0,ysol(xi)
do 10 i=2,nstep
x=xi+dx*float(i)
x1=x-dx
term1= 3.*x1*(y1-y0)/dx
term2=5.*y1
c print*,'term1,term2=',term1,term2
y2=2.*y1-y0+dx**2*(term1 +term2)/x1**2
if(i.eq.kount)then
print 100,x,y2,ysol(x)
kount=kount + kp
endif
y0=y1
y1=y2
10 continue
100 format(3(2x,e11.4))
stop
end
RUN
nstep,dx= 2000 0.000500000024
x ynum ysol
0.1000E+01 0.2000E+01 0.2000E+01
0.1025E+01 0.2107E+01 0.2107E+01
0.1050E+01 0.2228E+01 0.2229E+01
0.1075E+01 0.2365E+01 0.2366E+01
0.1100E+01 0.2519E+01 0.2520E+01
0.1125E+01 0.2690E+01 0.2691E+01
0.1150E+01 0.2880E+01 0.2881E+01
0.1175E+01 0.3089E+01 0.3091E+01
0.1200E+01 0.3319E+01 0.3322E+01
0.1225E+01 0.3572E+01 0.3575E+01
0.1250E+01 0.3849E+01 0.3852E+01
0.1275E+01 0.4150E+01 0.4154E+01
0.1300E+01 0.4478E+01 0.4482E+01
0.1325E+01 0.4834E+01 0.4839E+01
0.1350E+01 0.5219E+01 0.5225E+01
0.1375E+01 0.5636E+01 0.5642E+01
0.1400E+01 0.6085E+01 0.6093E+01
0.1425E+01 0.6569E+01 0.6578E+01
0.1450E+01 0.7090E+01 0.7099E+01
0.1475E+01 0.7649E+01 0.7660E+01
0.1500E+01 0.8249E+01 0.8260E+01
0.1525E+01 0.8891E+01 0.8904E+01
0.1550E+01 0.9578E+01 0.9592E+01
0.1575E+01 0.1031E+02 0.1033E+02
0.1600E+01 0.1109E+02 0.1111E+02
0.1625E+01 0.1193E+02 0.1195E+02
0.1650E+01 0.1282E+02 0.1284E+02
0.1675E+01 0.1376E+02 0.1378E+02
0.1700E+01 0.1476E+02 0.1479E+02
0.1725E+01 0.1583E+02 0.1585E+02
0.1750E+01 0.1696E+02 0.1698E+02
0.1775E+01 0.1815E+02 0.1818E+02
0.1800E+01 0.1942E+02 0.1945E+02
0.1825E+01 0.2076E+02 0.2079E+02
0.1850E+01 0.2217E+02 0.2221E+02
0.1875E+01 0.2367E+02 0.2371E+02
0.1900E+01 0.2525E+02 0.2529E+02
0.1925E+01 0.2691E+02 0.2695E+02
0.1950E+01 0.2866E+02 0.2871E+02
0.1975E+01 0.3050E+02 0.3056E+02
0.2000E+01 0.3244E+02 0.3250E+02