keepfloat : true$
load("stringproc")$
exact_soln_y(x) := block([],(2.0*x+1.0)^sin(x))$
solvediffeq(f,y,xl,x0,y0,n,xf,h) := block([yy,xx,true_y,t0,tf,est_time,iter,adj,series,series2],
  x : first(xl),
  iter : 0,
  xx : x0,
  yy : y0,
  t0 : elapsed_real_time(),
  series : taylor(first(f),first(xl),x0,n),
  series2 : integrate(series,x),
  adj : yy - subst(xx,x,series2),
  while (xx < xf and iter < 100) do (
    true_y : exact_soln_y(xx),
    print(iter," Iterations x = ",xx," y = ",yy," true y = ", true_y), 
    xx : xx + h,
    series : taylor(first(f),first(xl),xx,n),
    series2 : integrate(series,x),
    yy : subst(xx,x,series2) + adj, 
    true_y : exact_soln_y(xx),
    iter : iter + 1
    ),
  tf : elapsed_real_time(),
  if (xx - x0 > 0.0) then (
    est_time : (tf - t0) * (xf - x0) / (xx - x0)/60.0/60.0,
    printf(true,"Estimated Time = ~f Hours ~%",est_time)
    )
  else (
    printf(true,"Zero time")
    )
  );
solvediffeq([ ( 2.0 * x + 1.0) ^( sin ( x ) ) * ( cos ( x ) * log ( 2.0* x + 1.0 )+ ( 2.0 * sin ( x ) ) / ( 2.0 * x + 1.0 ))],[y] ,[x],0.1,exact_soln_y(0.1),30,1.0,0.00001)$