using Plots

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=100    # final time

function trapsolve(y0,t0,T,N)
	local Ys=zeros(N)
	local Ts=zeros(N)
	global yn,tn,h
	h=(T-t0)/N  # size of step
	yn=y0
	tn=t0
	for n=1:N
		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))
		# euler step
#		yn=yn+h*f(tn,yn)
		tn=t0+n*h
		Ys[n]=yn
		Ts[n]=tn
	end
	return Ts,Ys
end

Ts,Ys=trapsolve(y0,t0,T,64)
plot(Ts,Ys)