yexact(t)=1/(1/y0+sin(t0)-sin(t))
f(t,y)=y^2*cos(t)
dfdy(t,y)=2*y*cos(t)
phi(z)=yn+h/2*(f(tn,yn)+f(tn+h,z))-z
dphi(z)=h/2*dfdy(tn+h,z)-1
g(z)=z-phi(z)/dphi(z)

y0=0.8   # initial condition
t0=0     # initial time
T=2      # final time
N=32     # number of steps
h=(T-t0)/N  # size of step

yn=y0
tn=t0
for n=1:N
	global yn,tn
	local zn
	# Newton's method
	zn=yn
	for i=1:3
		zn=g(zn)
	end
	# trapezoid step
	yn=yn+h/2*(f(tn,yn)+f(tn+h,zn))
	tn=t0+n*h
end