# File diffeq.rb, line 341
def radius_real(y_higher,order_diff)
#local  m,n,rm0,rm1,hdrc,rcs,ord_no,ans
# Computes radius of convergence && order of pole from 3 adjacent Taylorv series terms.
# Applies to pole of arbitrary order on the real axis,
# Due to Prof. George Corliss.

ans = Array.new(3)
n = $max_taylor
m = n
while (m > 3) && ((y_higher[1][m].abs < $small_float) || (y_higher[1][m-1].abs < $small_float) || (y_higher[1][m-2].abs < $small_float )) do
if DEBUG == 1
$stderr.puts "radius_real m == " + m.to_s
$stderr.puts "m == " + m.to_s + "y[1][m] == " + (y_higher[1][m]).to_s +  " y[1][m-1] == " + (y_higher[1][m-1]).to_s + "y[1][m-2] == " + (y_higher[1][m-2]).to_s
end
m -= 1
end
if (m > 3) then
if DEBUG == 1 then
$stderr.puts "m == " + m.to_s + "y[1][m] == " + (y_higher[1][m]).to_s + "y[1][m-1] == " + (y_higher[1][m-1]).to_s + "y[1][m-2] == " + (y_higher[1][m-2]).to_s
end

rm0 = y_higher[1][m]/y_higher[1][m-1]
rm1 = y_higher[1][m-1]/y_higher[1][m-2]
hdrc = ap_int(m-1)*rm0-ap_int(m-2)*rm1
if DEBUG == 1
$stderr.puts "re hdrc == " + hdrc.to_s
end
if hdrc.not_zero($small_float) then
rcs = $h/hdrc
ord_no = ap_int(m-1)*rm0/hdrc -ap_int(m) + ap_int(2)
if DEBUG == 1
$stderr.puts "re rad == " + rcs.to_s
$stderr.puts "re order == " + ord_no.to_s
end
ans[1] = rcs
ans[2] = ord_no
return(ans)
else
ans[1] = $large_float
ans[2] = $large_float
return(ans)
end
else
ans[1] = $large_float
ans[2] = $large_float
return(ans)
end
return(ans)
end



def radius_complex(y_higher,order_diff)
#local m,rm0,rm1,rm2,rm3,rm4,rcs,ord_no,nr1,nr2,dr1,dr2,ds1,ds2,ans,REc,IMc,REordno,IMordno,n,cnt,rad_c
##
# computes radius of convergence for complex conjugate pair of poles.
# from 6 adjacent Taylor series terms
# Also computes order of poles.
# Due to Manuel Prieto.
# With a correction by Dennis J. Darland
ans = Array.new(3)
#DJDSTOP
#   if iter == 0 then
#     ans[1] = ($large_float)
#     ans[2] = ($large_float)
#     return(ans)
#   end
# #DJDSTART
n = $max_taylor
cnt = 0
while (cnt < 5  && n >= 1) do
if ((y_higher[1][n].abs) > ($small_float)) then
cnt = cnt + 1
else
cnt = 0
end

n -= 1
end
m = n + cnt
if (m < 6) then
ans[1] = ($large_float)
ans[2] = ($large_float)
return(ans)
elsif y_higher[1][m].abs >= ($large_float) || y_higher[1][m-1].abs >=($large_float) || y_higher[1][m-2].abs >= ($large_float) || y_higher[1][m-3].abs >= ($large_float) || y_higher[1][m-4].abs >= ($large_float) || y_higher[1][m-5].abs >= ($large_float) then
ans[1] = ($large_float)
ans[2] = ($large_float)
return(ans)
else
rm0 = (y_higher[1][m])/(y_higher[1][m-1])
rm1 = (y_higher[1][m-1])/(y_higher[1][m-2])
rm2 = (y_higher[1][m-2])/(y_higher[1][m-3])
rm3 = (y_higher[1][m-3])/(y_higher[1][m-4])
rm4 = (y_higher[1][m-4])/(y_higher[1][m-5])
if DEBUG == 1
$stderr.puts "rm0 == " + rm0.to_s
$stderr.puts "rm1 == " + rm1.to_s
$stderr.puts "rm2 == " + rm2.to_s
$stderr.puts "rm3 == " + rm3.to_s
$stderr.puts "rm4 == " + rm4.to_s
end
nr1 = ap_int(m-1)*rm0 - $two*ap_int(m-2)*rm1 + ap_int(m-3)*rm2
nr2 = ap_int(m-2)*rm1 - $two*ap_int(m-3)*rm2 + ap_int(m-4)*rm3
if DEBUG == 1
$stderr.puts "nr1 == " + nr1.to_s
$stderr.puts "nr2 == " + nr2.to_s
end
dr1 = $minus_one/rm1 + $two/rm2 - $one/rm3
dr2 = $minus_one/rm2 + $two/rm3 - $one/rm4
ds1 = $three/rm1 - $eight/rm2 + $five/rm3
ds2 = $three/rm2 - $eight/rm3 + $five/rm4
if DEBUG == 1
$stderr.puts "dr1 ==" +dr1.to_s
$stderr.puts "dr2 ==" + dr2.to_s
$stderr.puts "ds1 ==" + ds1.to_s
$stderr.puts "ds2 ==" + ds2.to_s
end
if (nr1 * dr2 - nr2 * dr1).abs <= $small_float || (dr1.abs) <= $small_float then
ans[1] = ($large_float)
ans[2] = ($large_float)
return(ans)
else

#term + dr1*dr2 added to Manuel Prieto's analysis
if (nr1*dr2 - nr2 * dr1).not_zero($small_float) then
rcs = ((ds1*dr2 - ds2*dr1 +dr1*dr2)/(nr1*dr2 - nr2 * dr1))
#              rcs = (ds1*dr2 - ds2*dr1)/(nr1*dr2 - nr2 * dr1)
ord_no = (rcs*nr1 - ds1)/($two*dr1) -ap_int(m)/$two
if rcs.not_zero($small_float) then # && rcs < $thousand * $h then
if DEBUG == 1
$stderr.puts "rcs == " + rcs.to_s
end
if rcs > $zero then
rad_c =( rcs.sqrt) * $h
else
rad_c =$large_float
end
else
rad_c = $large_float
ord_no = $large_float
end
else
rad_c = $large_float
ord_no = $large_float
end
$stderr.puts "radius == " + rad_c.to_s
$stderr.puts "order == " + ord_no.to_s
ans[1] = rad_c
ans[2] = ord_no
return(ans)
end
ans[1] = ($large_float)
ans[2] = ($large_float)
return(ans)
end
ans[1] = ($large_float)
ans[2] = ($large_float)
return(ans)

end
# ans[1] = ($large_float)
# ans[2] = ($large_float)
# return(ans)
# end

def which_radius(prev_pole,real_pole,complex_pole,print_flag)
#local pole
pole = Array.new(3)
$point_9 = Apfp.new(9,-1,50,-NUM_DIGITS)
if (real_pole[1] == ($large_float) || real_pole[2] == ($large_float)) && (complex_pole[1] != ($large_float)) && (complex_pole[2] != ($large_float)) && (complex_pole[1] > $point_9) && (complex_pole[2] > ($zero))  then
pole[1] = complex_pole[1]
pole[2] = complex_pole[2]
if $correct_start_flag != 1 && print_flag == 1 then
puts "complex estimate of poles used"
end
if $correct_start_flag != 1 && print_flag == 1 then
$stderr.puts "complex estimate of poles used"
end
elsif (real_pole[1] != $large_float) && (real_pole[2] != $large_float) && (real_pole[1] > $point_9) && (real_pole[2] > $zero) && ((complex_pole[1] == $large_float) || (complex_pole[2] == $large_float) || (complex_pole[1] <= $point_9 ) || (complex_pole[2] <= $zero))
pole[1] = real_pole[1]
pole[2] = real_pole[2]
if $correct_start_flag != 1 && print_flag == 1 then
puts "real estimate of pole used"
end
if $correct_start_flag != 1 && print_flag == 1 then
$stderr.puts "real estimate of pole used"
end
elsif ((real_pole[1] == $large_float) || (real_pole[2] == $large_float)) && ((complex_pole[1] == $large_float) || (complex_pole[2] == $large_float))
pole[1] = $large_float
pole[2] = $large_float
if $correct_start_flag != 1 && print_flag == 1 then
puts "NO POLE 1"
end
if $correct_start_flag != 1 && print_flag == 1 then
$stderr.puts "NO POLE 1"
end
elsif (real_pole[1] < complex_pole[1]) && (real_pole[1] > $point_9) && (real_pole[2] > $zero)
pole[1] = real_pole[1]
pole[2] = real_pole[2]
if $correct_start_flag != 1 && print_flag == 1 then
puts "real estimate of pole used"
end
if $correct_start_flag != 1 && print_flag == 1 then
$stderr.puts "real estimate of pole used"
end
elsif (complex_pole[1] != $large_float) && (complex_pole[2] != $large_float) && (complex_pole[1] > $point_9) && (complex_pole[2] > $zero)
pole[1] = complex_pole[1]
pole[2] = complex_pole[2]
if $correct_start_flag != 1 && print_flag == 1 then
puts "complex estimate of poles used"
end
if $correct_start_flag != 1 && print_flag == 1 then
$stderr.puts "complex estimate of poles used"
end
else
pole[1] = $large_float
pole[2] = $large_float
if $correct_start_flag != 1 && print_flag == 1 then
puts "NO POLE 2"
end
if $correct_start_flag != 1 && print_flag == 1 then
$stderr.puts "NO POLE 2"
end
end
if (pole[1] != $large_float) && ($correct_start_flag != 1 && print_flag == 1) then
puts "radius of convergence == " + pole[1].to_s
puts "order of singularity  == " + pole[2].to_s
$stderr.puts "radius of convergence == " + pole[1].to_s
$stderr.puts "order of singularity  == " + pole[2].to_s
end
return(pole)
end


def compute_order(log10norm, log10abserr, log10relerr)
#*
# * this is to determine the 'optimal' degree of the taylor expansion.
# *
# * parameters:
# * log10norm: base 10 log of the norm of the initial condition
# * log10abserr: base-10 log of the absolute error required
# * log10relerr: base-10 log of the relative error required
# *
# * returned value: 'optimal' degree.
#*/
# When I studied it I couldn't really understand Angel Jorba && Maorong Zou 's
# method for determining the order.  I don't adjust it anymore.
# It will be set to the problem file max_terms which defaults to 30 which seems
# a good value most of the time.
# I've retained $max_taylor in case I find a better way.
return($max_taylor)
end

def compute_step_size(x,print_flag,log10norm,log10abserr,log10relerr)
#/*
# * it looks for a step size to keep the truncation error (both
# * absolute and relative) below the absolute error tol, and also tries
# * to reduce cancellations of big terms in the summation of the taylor
# * series.
# *
# * WARNING: there are two versions of this function. this is the one
# * that uses MY_FLOAT variables for all the computations
# * */

#local dh,hd,j,aaa,r2,cc,hd2,i
#local  h,r,z,hs
cs_info("begin compute step size")
if DEBUG == 1 then
puts "compute_step_size last time = " + $adjusted_h_last_time.to_s + "optimal = " + $reached_optimal_h.to_s + "start_flag = " + $correct_start_flag.to_s
$stderr.puts "compute_step_size last time = " + $adjusted_h_last_time.to_s + "optimal = " + $reached_optimal_h.to_s + "start_flag = " + $correct_start_flag.to_s
end
if $adjusted_h_last_time == 1 || $reached_optimal_h == 1
$adjusted_h_last_time = 0
return $h
end
#  we will use the array norms to store the norms of the normalized
#  derivatives
#*/
#/*
#  we compute the step size according to the first algorithm
#*/
cs_info("begin compute step size A")
get_norms()
dh=compute_step_size_A(x,print_flag,log10norm,log10abserr,log10relerr)
cs_info("After compute step size A")

h = (dh)
hs = (h)
if DEBUG == 1
$stderr.puts "dh == h == hs ==" + dh.to_s
end
#/*
#  next loop runs over the norms of the taylor coefficients, checking
#  whether the summation of the taylor series (from low to high order)
#  involves increasing terms. in this case, it reduces the stepsize
#  accordingly.
#*/
hd = h
z = $norms[1]*hd
hdj = hd
flag = 0
j = 1
while j<= $max_taylor-1 do
hdj *= hd
r= $norms[j+1]*hdj
if DEBUG == 1
$stderr.puts "loop over norms j == " + j.to_s + "norm == " + r.to_s
end
if (r <= z) then
j += 1
next
end
if DEBUG == 1
$stderr.puts "r == " + r.to_s
$stderr.puts "z == " + z.to_s
end
# reduce hd (& hdj)
flag = 1
hdj = z/$norms[j+1]
a = (r/z) ** j
if DEBUG == 1
$stderr.puts "a == " + a.to_s
end
if DEBUG == 1
$stderr.puts "cc == " + cc.to_s
$stderr.puts "hd == " + hd.to_s
end
hd2 = hd
hd=hd2/a
#    z = (r)
if DEBUG == 1
$stderr.puts "order " + j.to_s + " reducing h from " + hd2.to_s + " to " + hd.to_s
end
j += 1
end
if DEBUG == 1
$stderr.puts "compute_step_size h == " + h.to_s
$stderr.puts "compute_step_size hd == " + hd.to_s
end
if hd < ($hmin_init) then
return c0($hmin_init)
end
if hd > ($hmax) then
return c0($hmax)
end
return c0(hd)
end

def estimate_trunc_err
get_norms()
d1 = $norms[$max_taylor-$max_order]
d2 = $norms[$max_taylor-1-$max_order]
d3 = $norms[$max_taylor-2-$max_order]
d4 = $norms[$max_taylor-3-$max_order]

return (d1.abs+d2.abs+d3.abs+d4.abs)
end


def compute_step_size_A(x,print_flag,log10norm,log10abserr,log10relerr)

#local sz,d1,d2,adj,sz2,i,d_min,log10d1,log10d2,r,nt,log10ro,log10ro1,log10ro2,found,flag1,flag2,log10err,f1,f2,log10f1,log10f2

if log10abserr < (log10relerr + log10norm) then
log10abserr = log10relerr + log10norm
end
d1 = ap_max($norms[$max_taylor-$max_order],$small_float)
d2 = ap_max($norms[$max_taylor-1-$max_order],$small_float)
d3 = ap_max($norms[$max_taylor-2-$max_order],$small_float)
d4 = ap_max($norms[$max_taylor-3-$max_order],$small_float)
if DEBUG == 1
$stderr.puts "d1 == " + d1.to_s
$stderr.puts "d2 == " + d2.to_s
$stderr.puts "d3 == " + d3.to_s
$stderr.puts "d4 == " + d4.to_s
$stderr.puts "log10abserr == " + log10abserr.to_s
end
log10d1 = d1.about_log10
log10d2 = d2.about_log10
log10d3 = d3.about_log10
log10d4 = d4.about_log10
r = log10abserr
log10ro4 = (r - log10d4)/ap_int($max_taylor-4-$max_order)
log10ro3 = (r - log10d3)/ap_int($max_taylor-3-$max_order)
log10ro2 = (r - log10d2)/ap_int($max_taylor-2-$max_order)
log10ro1 = (r - log10d1)/ap_int($max_taylor-1-$max_order)
log10roa = ap_max(log10ro1,log10ro2)
log10rob = ap_max(log10ro3,log10ro4)
log10ro = ap_max(log10roa,log10rob)
if DEBUG == 1 then
$stderr.puts "r ==" + r.to_s
$stderr.puts "log10d1==" + log10d1.to_s
$stderr.puts "log10d2==" + log10d2.to_s
$stderr.puts "log10d3==" + log10d3.to_s
$stderr.puts "log10ro1==" + log10ro1.to_s
$stderr.puts "log10ro2==" + log10ro2.to_s
$stderr.puts "log10ro3==" + log10ro3.to_s
$stderr.puts "log10ro==" + log10ro.to_s
end
#  sz = fix_h($ten ** (log10ro - $four) * $h)
#  sz = fix_h(Apfp.new(10 ** ((log10ro.mant * (10 ** log10ro.expt) - ($max_order * $no_eqs)),0,50,-NUM_DIGITS).norm * $h)
newmant = 10 ** ((log10ro.mant * (10 ** log10ro.expt)) - (($tot_order + 1) * $no_eqs) + 1)
#  newmant = 10 ** (log10ro.mant * (10 ** log10ro.expt))
sz = (Apfp.new(newmant,0,50,-NUM_DIGITS+newmant.to_s.size).norm * $h).norm
sz = fix_h(sz.clone)
if DEBUG == 1 then
$stderr.puts "sz(computed new $h) == " + sz.to_s
end
cs_info('After sz computed')
if $hmax < sz && $h != $hmin_init
$unchanged_h_cnt = 0
$reached_optimal_h = 1
$optimal_clock_start_sec = time_sec
$optimal_start  = x[1]
if $hmax < sz then
sz= $hmax
end
end
if $h * $two < sz then
puts "$h suggested at " + ($h * $two).to_s
$stderr.puts "$h suggested at " + ($h * $two).to_s
puts "Requested accuracy was $log10abserr == " + $log10abserr.to_s
$stderr.puts "Requested accuracy was $log10abserr == " + $log10abserr.to_s
cs_info("$h suggested ")
return c0($two * $h)
end
return(sz)
end


def adjust_step(x,y_higher,hold,hnew,order_diff)
#local order,term_no,accum,r
if DEBUG == 1
$stderr.puts "in adjust hold == " + hold.to_s
$stderr.puts "in adjust hnew == " + hnew.to_s
end
if hnew < $hmin_init then
hnew = $hmin_init
end
if hold >= hnew * $two then
puts "Adjusting $h from " + hold.to_s + " to " + hnew.to_s
$stderr.puts "Adjusting $h from " + hold.to_s + " to " + hnew.to_s
puts "Requested accuracy was $log10abserr == " + $log10abserr.to_s
$stderr.puts "Requested accuracy was $log10abserr == " + $log10abserr.to_s
#$x[2]= hnew
accum= $one
r = hnew/hold
if DEBUG == 1
$stderr.puts "r == " + r.to_s
end

term_no = 1
while term_no <= $max_terms do
order = 1
while order <= order_diff+1 do
if DEBUG == 1
$stderr.puts "accum == " + accum.to_s
end
y_higher[order][term_no] = y_higher[order][term_no] * accum
if DEBUG == 1
$stderr.puts "y_higher[" + order.to_s + "][" + term_no.to_s + "] == " + (y_higher[order][term_no]).to_s
end
order += 1
end
accum = accum*r
term_no += 1
end
return (c0(hnew))
else
return c0(hold)
end
end
def factorial_1(n)
# $stderr.puts "factorial " + n.to_s
return ap_int($fact_tbl[n])
end
def array(i,j)
#local L
l = Array.new(i)
k = 1
while k <= j do
l[k] = Array.new(j)
k += 1
end
return l
end

def Si(x)
puts "Si not implemented in Icon"
exit(0)
end
def erf(x)
puts "erf not implemented in Icon"
exit(0)
end

def time_sec()
t = Time.now
return t.to_f
end

def wrk_prod(f,n)
#local p,i,r
p = f
r = f - 1
i = n
while (i >= 1 && r >= 1) do
p *= r
r -= 1
i -= 1
end
return ap_int(p)
end

def chk_data
#local errflag
errflag = 0

if $max_terms < 30 || $max_terms > 512 then
puts "Illegal max_terms: " + $max_terms.to_s + " -- Using 30"
$stderr.puts "Illegal max_terms: " + $max_terms.to_s
$max_terms = 30
$max_taylor = 30
end
if ($look_poles != 0 && $look_poles != 1) then
puts "Illegal look_poles"
$stderr.puts "Illegal look_poles"
errflag = 1
end
if ($adjust_h != 0 && $adjust_h != 1) then
puts "Illegal adjust_h "
$stderr.puts "Illegal adjust_h"
errflag = 1
end
if ($look_poles == 1 && $adjust_h == 0) then
puts "Illegal look_poles or adjust_h"
$stderr.puts "Illegal look_poles or adjust_h"
errflag = 1
end
if ($max_iter < 10) then
puts "Illegal max_iter"
$stderr.puts "Illegal max_iter"
errflag = 1
end
if errflag == 1 then
exit
end
end
def comp_expect(t_end2,t_start2,t2,clock_sec)
#local hours,ms,ms2,d,sub1,sub2
ms2 = clock_sec
d = 3600.0
sub1 = (t_end2-t_start2).to_f
sub2 = (t2-t_start2).to_f
if sub2 != 0.0
hours = sub1*ms2/(sub2*d)
else
hours = 0.0
end
return hours
end
def comp_expect_sec(t_end2,t_start2,t2,clock_sec)
#local hours,ms,ms2,d,sub1,sub2
ms2 = clock_sec
sub1 = (t_end2-t_start2)
sub2 = (t2-t_start2)
if sub2.not_zero($small_float) then
r = (sub1/sub2).to_f
sec_left = r * ms2
else
sec_left = 0.0
end

return sec_left
end
def time_elapsed(clock_sec_v)
#local ms,hrs
sec = 1.0*clock_sec_v
hrs = sec/(3600.0)
return hrs.to_f
end
def fix_h(h)
mant = h.mant.to_s
expt = h.expt
expt = expt.to_i + mant.to_s.size - 2
mant = mant[0..1]
errmant = 5
errexpt = -NUM_DIGITS+expt+2
return Apfp.new(mant.to_i,expt.to_i,errmant.to_i,errexpt.to_i)
end
def sec_to_display(sec)
if sec == 0.0 then
return " Unknown"
end

sec_in_millinium =  SEC_IN_MIN *  MIN_IN_HOUR * HOURS_IN_DAY * DAYS_IN_YEAR * YEARS_IN_CENTURY * CENTURIES_IN_MILLINIUM
milliniums = sec/sec_in_millinium
millinium_int = milliniums.to_i
centuries = (milliniums - millinium_int.to_f)* CENTURIES_IN_MILLINIUM
cent_int = centuries.to_i
years = (centuries-cent_int.to_f) * YEARS_IN_CENTURY
years_int = years.to_i
days = (years - years_int.to_f)  * DAYS_IN_YEAR
days_int = days.to_i
hours = (days - days_int.to_f) * HOURS_IN_DAY
hours_int = hours.to_i
minutes = (hours - hours_int.to_f) * MIN_IN_HOUR
minutes_int = minutes.to_i
seconds = (minutes - minutes_int.to_f) * SEC_IN_MIN
sec_int = seconds.to_i
str = " "
if millinium_int > 0 then
str = str + millinium_int.to_s + " Millinia "
end
if cent_int > 0 then
str = str + cent_int.to_s + " Centuries "
end
if years_int > 0 then
str = str + years_int.to_s + " Years "
end
if days_int > 0 then
str = str + days_int.to_s + " Days "
end
if hours_int > 0 then
str = str + hours_int.to_s + " Hours "
end
if minutes_int > 0 then
str = str + minutes_int.to_s + " Minutes "
end
if sec_int > 0 then
str = str + sec_int.to_s + " Seconds "
end
return str
end


def cs_info(str)
if DEBUG == 1
puts "cs_info " + str + " $correct_start_flag == " + $correct_start_flag.to_s + "$h = " + $h.to_s + "$reached_optimal_h = " + $reached_optimal_h.to_s
$stderr.puts "cs_info " + str + " $correct_start_flag == " + $correct_start_flag.to_s + "$h = " + $h.to_s + "$reached_optimal_h = " + $reached_optimal_h.to_s
end
end

# END ATS LIBARY BLOCK
# BEGIN USER DEF BLOCK
def exact_soln_y(y_str,x)
case y_str
when 'y'
return($minus_one*(x.cos))
else
return $zero
end

end

# END USER DEF BLOCK
# END OUTFILE5
# BEGIN OUTFILEMAIN
 $orig_start_sec = time_sec
$curr_iter_when_opt = 0
rconst = $ApConst.new
rconst.one.setrconst(rconst)
rconst.init2
rconst.one.setrconst(rconst)
$RConst = rconst
cconst = ApcConst.new(rconst)
cconst.one.setcconst(cconst)
cconst.one.setrconst(rconst)
$CConst = cconst
set_constants
 $max_terms = 30
 $max_taylor = 30
 print_flag = 1
 $min_terms = 30
 $max_order = 1
 $no_eqs = 1
 $tot_order = 1
# BEGIN CONST BLOCK
$max_terms = 30
# END CONST BLOCK
iter = -1
$norms = Array.new($max_terms + 1)
$max_iter = 50000
$max_hours = 0.0
$max_minutes = 15.0
$tmp0 = Array.new($max_terms + 1)
set_z($tmp0)
$tmp1 = Array.new($max_terms + 1)
set_z($tmp1)
$tmp2 = Array.new($max_terms + 1)
set_z($tmp2)
$x = Array.new($max_terms + 1)
set_z($x)
$tmp1_g = Array.new($max_terms + 1)
set_z($tmp1_g)
$y = Array.new($max_terms + 1)
set_z($y)
$y_init = Array.new($max_terms + 1)
$y_higher = array(2,$max_terms+1)
$y_higher_work = array(2,$max_terms+1)
puts "##############ECHO OF PROBLEM#################"
$stderr.puts "##############ECHO OF PROBLEM#################"
puts "##############" + "sin.ode" +"#################"
$stderr.puts "###############" + "sin.ode" + "}#################"
puts "diff ( y , x , 1 ) =  sin ( x ) ;"
$stderr.puts "diff ( y , x , 1 ) =  sin ( x ) ;"
puts "!"
$stderr.puts "!"
puts "$max_terms = 30"
$stderr.puts "$max_terms = 30"
puts "!"
$stderr.puts "!"
puts "$x_start = 0.0"
$stderr.puts "$x_start = 0.0"
puts "$x_end = 10.0 "
$stderr.puts "$x_end = 10.0 "
puts "$y_init[1] = exact_soln_y('y',c($x_start))"
$stderr.puts "$y_init[1] = exact_soln_y('y',c($x_start))"
puts "$hmax = 0.01 "
$stderr.puts "$hmax = 0.01 "
puts "$adjust_h = 1"
$stderr.puts "$adjust_h = 1"
puts "$look_poles = 0"
$stderr.puts "$look_poles = 0"
puts "$max_iter = 200000"
$stderr.puts "$max_iter = 200000"
puts "$log10relerr = -8"
$stderr.puts "$log10relerr = -8"
puts "$log10abserr = -8"
$stderr.puts "$log10abserr = -8"
puts "!"
$stderr.puts "!"
puts "def exact_soln_y(y_str,x)"
$stderr.puts "def exact_soln_y(y_str,x)"
puts "case y_str"
$stderr.puts "case y_str"
puts "  when 'y'"
$stderr.puts "  when 'y'"
puts "          return($minus_one*(x.cos))"
$stderr.puts "          return($minus_one*(x.cos))"
puts "  else"
$stderr.puts "  else"
puts "          return $zero"
$stderr.puts "          return $zero"
puts "  end"
$stderr.puts "  end"
puts ""
$stderr.puts ""
puts "end"
$stderr.puts "end"
puts ""
$stderr.puts ""
puts "!"
$stderr.puts "!"
puts "#######END OF ECHO OF PROBLEM#################"
$stderr.puts "#######END OF ECHO OF PROBLEM#################"
set_z($norms)
set_z($y_init) 
set_zz($y_higher,1)
set_zz($y_higher_work,1)
$correct_start_flag = 1
$unchanged_h_cnt = 0
$adjust_h = 1
$look_poles = 0
$warned = 0
$warned2 = 0
$small_float = c0(1.0e-200)
$smallish_float = c0(1.0e-10)
$large_float = c0(1.0e100)
$almost_1 = c0(0.99)
$progress_ind = 1
$log10abserr = -8
$log10relerr = -8
$hmax = c0(0.01)
$old_log10eps=c0(-1.0)
$old_nt=0
$adjusted_h_last_time = 0
$