#Instructions: use this function to compute approximate numerical values of 
#solutions to initial value problems using the backward euler method
#
#if you want to solve dy/dt=f(t,y) for y(t0)=y0 with step size h,
#try backeulerapprox(f,t0,y0,h,n) where n is the number of values you want to see
#
#for example, if y'=3+t-y and y(0)=1, and step size h=0.1, try
#backeulerapprox(3+t-y,0,1,0.1,4);
#
backeulerapprox:=proc(f,t0,y0,h,n) local yprev, ycurr, tprev, tcurr, i:
ycurr:=y0:
tcurr:=t0:

for i from 1 to n do
yprev:=ycurr:
tprev:=tcurr:

#If y_n=yprev and t_n=tprev, then...

#y_(n+1)=y_n +f(t_(n+1),t_(y+1))*h 
#= y_n+f(t_n+h,y_(n+1)*h
#notice that this is an equation with y_(n+1) on BOTH sides,
#so we must SOLVE it for y_(n+1)
#(this is equation 13 on page 445)
ycurr:=solve(y=yprev+subs({t=tprev+h},f)*h,y):

#t_(n+1)=t_n+h
tcurr:=tprev+h:

print("when t=", tcurr, ", y is approximately =", ycurr):
od:

end: