subroutine pcstep(x0,x1,h,y0,y1,y,n,f,err,f1,f2,c,yt)
dimension y0(n),y1(n),y(n),c(n),yt(n),f1(n),f2(n)
external f

x2=x1+h; call f(x1,y1,f1)
do i=1,n
	c(i)=(y0(i)-y1(i)-f1(i)*(x0-x1))/(x0-x1)**2
	yt(i)=y1(i)+f1(i)*h+c(i)*h*h; end do
call f(x2,yt,f2)
err=0
do i=1,n
	c(i)=(f2(i)-f1(i)-2*c(i)*(x2-x1))/(2*(x2-x1)*(x2-x0)+(x2-x1)**2)
	corr=c(i)*(x2-x1)**2*(x2-x0)
	err=max(err,abs(corr))
	y(i)=yt(i)+corr;
end do
err=err/2
return; end