import sys
from sage.all import *
import cProfile
import random
import os
import math
import numpy
import time
import pprint



def main():
    problem_1(skip1k=True) # otherwise takes a few hours
    problem_2()
    problem_3()
    problem_4()


###

def get_dydt(y,t):
    to_ret = -y*t + 0.1*y**3
    return to_ret


def get_soln(y0, t0, tmax, deltat, F=None, spacing=1, with_deriv=False):
    y_vals = {t0:y0}
    yderiv_vals = {}
    y = float(y0)
    t = float(t0)
    counter = 0
    while t < tmax:
        if F is None:
            dydt = get_dydt(y,t)
        else:
            dydt = F(y,t)
        if with_deriv and ((counter%spacing) == 0):
            yderiv_vals[t] = dydt
        new_t = t + deltat
        new_y = y + dydt * deltat
        t = new_t
        y = new_y
        if (counter%spacing) == 0:
            y_vals[t] = y
        counter += 1
    y_vals[t] = y
    if with_deriv:
        yderiv_vals[t] = dydt
        to_ret = [y_vals, yderiv_vals]
    else:
        to_ret = y_vals
    return to_ret


def pts_plot_single(y_vals, spacing=1, plot_spacing=1, kmin=None, kmax=None, vmin=None, vmax=None, absvmax=None, with_deriv=False, **kwargs):
    def condition(k,v):
        kmin_cond = (kmin is None) or (k >= kmin)
        kmax_cond = (kmax is None) or (k <= kmax)
        vmin_cond = (vmin is None) or (v >= vmin)
        vmax_cond = (vmax is None) or (v <= vmax)
        absvmax_cond = (absvmax is None) or (abs(v) <= absvmax)
        to_ret = kmin_cond and kmax_cond and vmin_cond and vmax_cond and absvmax_cond
        return to_ret

    if type(y_vals) == tuple:
        y0, t0, tmax, deltat = y_vals
        if with_deriv:
            y_vals, yderiv_vals = get_soln(y0, t0, tmax, deltat, spacing=spacing, with_deriv=with_deriv)
        else:
            y_vals = get_soln(y0, t0, tmax, deltat, spacing=spacing)

    if 'size' not in kwargs:
        kwargs['size'] = 1
    
    pts = points([(k,v) for k,v in sorted(list(y_vals.items()))[::plot_spacing] if condition(k,v)], dpi=400, **kwargs) # sage only
    if with_deriv:
        pts_deriv = points([(k,v) for k,v in sorted(list(yderiv_vals.items()))[::plot_spacing] if condition(k,v)], dpi=400, **kwargs)
        pts = pts + pts_deriv
    return pts


def pts_plot_range(ymin, ymax, num_ys, t0, tmax, deltat, F=None, spacing=1, kmin=None, kmax=None, vmin=None, vmax=None, absvmax=None, **kwargs):
    deltay = float(ymax - ymin)/num_ys
    y0list = [ymin + deltay*i for i in range(num_ys+1)]
    soln_dict = {y0:get_soln(y0, t0, tmax, deltat, F=F, spacing=spacing) for y0 in y0list}
    pts_list = [pts_plot_single(soln_dict[y0], spacing=spacing, kmin=kmin, kmax=kmax, vmin=vmin, vmax=vmax, absvmax=absvmax, color=hue(3.0*i/len(soln_dict))) for i,y0 in enumerate(sorted(list(soln_dict.keys())))]
    return pts_list


###


def problem_1(skip1k=False):
    def F(z,x):
        return cos(100*x**2 + z**3)
    z0 = 1.0
    x0 = 0.0
    xmax = 10.0
    # To know how small deltax needs to be taken to get a reasonable value of z(10), I start by looking at how quickly the function F will be oscillating around x = 10.
    # We can think of the argument as (100x)x + z^3, and imagine that the angular frequency is on the order of 100x.
    # For x = 10, this angular frequency will be about 10^3,
    # so I definitely want deltax*10^3 to be significantly smaller than 1 (because otherwise I'll be "jumping to some random point on the cosine").
    # I have to worry about the z^3 term too. The spirit of this is to just try things, see what happens, and have that inform further reasoning,
    # so it's reasonable to take deltax = 10^(-5) and move on for now.
    #
    # Alternatively, a more principled approach is to consider a Taylor series expansion, written out to second order:
    # z(x + deltax) = z(x) + z'(x)*deltax + 1/2*z''(x)*deltax^2 + ...
    # Euler's method, and e.g. the code in get_soln, hinges on the linear approximation z(x + deltax) ~= z(x) + z'(x)*deltax
    # For this to be a reasonable approximation, we need z'(x)*deltax >> z''(x)*deltax^2. When will this be the case?
    # We have z'(x) = cos(100x^2 + z^3), and can compute z''(x) = -(200x + 3z^2*z')*sin(100x^2 + z^3).
    # We see that z'(x) is of size 1 (it's the cos of something), and z''(x) is of size about 200x + 3z^2.
    # It'll be clear later that, in our case, 0 < z < 1, so ultimately z'(x) ~= 1, and z''(x) ~= 200x
    # If we want z'(x)*deltax >> 1/2*z''(x)*deltax^2, then we need 1 >> 100x*deltax, i.e. deltax << 1/100x.
    # Going up to x=10, this means that deltax = 10^(-3) probably won't work, deltax = 10^(-4) is borderline, and deltax = 10^(-5) is likely fine.
    deltax = 10.0**(-5)

    ### a)
    # with deltax = 10^(-5) and xmax - x0 = 10^1, this will be about 10^6 steps, which should run in at most a couple seconds.
    start_time = time.time()
    z_vals = get_soln(z0, x0, xmax, deltax, F=F) # passing a function is pretty python-y, but see the function get_soln for an alternative (basically hardcoding it).
    end_time = time.time()
    runtime = end_time - start_time
    z10_val = sorted(list(z_vals.items()))[-1] # one of the disavantages of using a dict with float keys is that it's unwieldy to access specific values
    print('last (x,z): ' + str(z10_val))  # last (x,z): (10.00000999979059, 0.9440721176773775)
    print('runtime: ' + str(runtime) + ' seconds') # runtime: 0.8546097278594971 seconds
    # We should check that our answer appears to be reasonably accurate
    # One good starting point is to look at a plot of your numerical solution.
    # Before doing that, let's also try a deltax that won't be precise enough, so that we can compare.
    # Because the cosine's oscillations are probably something like 1000x around x = 10,
    # taking deltax = 0.001 means that the argument of the cosine will increase by an amount of order 1 each step, which is too big and will give random values.
    deltax_toobig = 10.0**(-3)
    start_time_bad = time.time()
    z_vals_bad = get_soln(z0, x0, xmax, deltax_toobig, F=F)
    end_time_bad = time.time()
    runtime_bad = end_time_bad - start_time_bad
    z10_val_bad = sorted(list(z_vals_bad.items()))[-1]
    print('last (x,z) for deltax_toobig: ' + str(z10_val_bad)) # last (x,z) for deltax_toobig: (10.000999999999896, 0.9453853178390224)
    print('runtime_bad: ' + str(runtime_bad) + ' seconds') # runtime_bad: 0.007567167282104492 seconds
    
    # the following is sage-only, sorry. Sagemath.org is where you can get this.
    plot_good = pts_plot_single(z_vals)
    plot_bad = pts_plot_single(z_vals_bad, color='red')
    plot_good_zoomed = pts_plot_single({x:z for x,z in z_vals.items() if x > 9.9})
    plot_bad_zoomed = pts_plot_single({x:z for x,z in z_vals_bad.items() if x > 9.9}, size=10, color='red')
    ga = graphics_array([plot_good, plot_bad, plot_good_zoomed, plot_bad_zoomed], 2, 2)
    ga.show(figsize=[22,16]) # a4f1.png
    # while the bad plot has nice patterns, it doesn't make sense as a continuous function: the function values are jumping all over the place.
    
    # Another way to get a sense of whether or not your function has converged is to compute it again with a different precision and see if you get the same answer
    deltax_moreprec = 10.0**(-6)
    start_time_moreprec = time.time()
    z_vals_moreprec = get_soln(z0, x0, xmax, deltax_moreprec, F=F)
    end_time_moreprec = time.time()
    runtime_moreprec = end_time_moreprec - start_time_moreprec
    z10_val_moreprec = sorted(list(z_vals_moreprec.items()))[-1]
    print('last (x,z) moreprec: ' + str(z10_val_moreprec)) # last (x,z) moreprec: (10.000000999267517, 0.9440720690646837)
    print('runtime_moreprec: ' + str(runtime_moreprec) + ' seconds') # runtime_moreprec: 10.820501804351807 seconds
    
    vals_to_check = sorted([(k,v) for k,v in z_vals.items() if k <= 10], reverse=True)[:20*10**2:10**2]
    vals_to_check_moreprec = sorted([(k,v) for k,v in z_vals_moreprec.items() if k <= 10], reverse=True)[:20*10**3:10**3]
    vals_to_check_bad = sorted([(k,v) for k,v in z_vals_bad.items() if k <= 10], reverse=True)[:20]
    print('deltax:   %.3E' % deltax + ' '*41 + '%.3E' % deltax_moreprec + ' '*51 + '%.3E' % deltax_toobig)
    for i in range(len(vals_to_check)):
        v1str = str(vals_to_check[i])
        v2str = str(vals_to_check_moreprec[i])
        v3str = str(vals_to_check_bad[i])
        print(v1str + ' '*(60 - len(v1str)) + v2str + ' '*(60 - len(v2str)) + v3str)
    # deltax:   1.000E-05                                         1.000E-06                                                   1.000E-03
    # (9.99999999979059, 0.9440761841122429)                      (9.999999999267517, 0.9440724757181846)                     (9.999999999999897, 0.9457877680671167)
    # (9.998999999790628, 0.9449094389456116)                     (9.998999999268266, 0.9449046048498061)                     (9.998999999999898, 0.946454108638366)
    # (9.997999999790666, 0.9446731699091457)                     (9.997999999269014, 0.9446756142745795)                     (9.997999999999898, 0.945498130225058)
    # (9.996999999790704, 0.9440380843523466)                     (9.996999999269763, 0.9440355986303942)                     (9.996999999999899, 0.9456282942441344)
    # (9.995999999790742, 0.9448043460141824)                     (9.995999999270511, 0.9447986847322744)                     (9.9959999999999, 0.9464772667167446)
    # (9.99499999979078, 0.9448046996948732)                      (9.99499999927126, 0.9448066042806008)                      (9.9949999999999, 0.9456427235705837)
    # (9.993999999790818, 0.9440397683763839)                     (9.993999999272008, 0.9440385639838582)                     (9.9939999999999, 0.9454883675929151)
    # (9.992999999790856, 0.9446765688726836)                     (9.992999999272756, 0.9446703832153793)                     (9.9929999999999, 0.9464520474744722)
    # (9.991999999790893, 0.9449150030932363)                     (9.991999999273505, 0.9449160564287822)                     (9.991999999999901, 0.9458080024593074)
    # (9.990999999790931, 0.9440820577960058)                     (9.990999999274253, 0.9440820893570003)                     (9.990999999999902, 0.9453799198530681)
    # (9.98999999979097, 0.9445363514973307)                      (9.989999999275001, 0.944529992709081)                      (9.989999999999903, 0.9463797623249508)
    # (9.988999999791007, 0.9449947026031381)                     (9.98899999927575, 0.9449946584323905)                      (9.988999999999903, 0.9459807630388373)
    # (9.987999999791045, 0.9441624077506493)                     (9.987999999276498, 0.9441635260297327)                     (9.987999999999904, 0.945312646150428)
    # (9.986999999791083, 0.9443953015043455)                     (9.986999999277247, 0.9443891416561018)                     (9.986999999999904, 0.9462657306095104)
    # (9.98599999979112, 0.9450367096501623)                      (9.985999999277995, 0.9450354096382384)                     (9.985999999999905, 0.946146719036919)
    # (9.984999999791158, 0.9442749663012712)                     (9.984999999278743, 0.9442769273219602)                     (9.984999999999905, 0.9452930927142741)
    # (9.983999999791196, 0.944265466527375)                      (9.983999999279492, 0.9442598679757919)                     (9.983999999999906, 0.9461189811864977)
    # (9.982999999791234, 0.9450369316765329)                     (9.98299999928024, 0.945034321387336)                       (9.982999999999906, 0.9462916912127006)
    # (9.981999999791272, 0.9444109278169068)                     (9.981999999280989, 0.9444134107222423)                     (9.981999999999907, 0.945323957416496)
    # (9.98099999979131, 0.9441583196527965)                      (9.980999999281737, 0.9441536039456326)                     (9.980999999999908, 0.9459516146858059)
    
    # We see that we get basically the same answer having used a deltax ten times smaller, so we can guess that the values we have are correct
    # We also see that the values for deltax = 0.001 don't match (look carefully!).
    # As we could have guessed from the plots and the second order Taylor series reasoning, deltax = 10^(-3) wasn't enough precision to be accurate.

    ### b)
    # From the plot, we immediately see a function oscillating more and more quickly, as we might guess from the differential equation saying its derivative oscillates quickly.
    # A notable feature is the downward trend, which we'll now speculate about:
    # The argument 100x^2 + z^3 will always be increasing in the range we're looking at, because 100x^2 is much bigger than our z, which is between 0 and 1.
    # Individually, 100x^2 is always increasing, but z^3 may be increasing or decreasing.
    # If z^3 is increasing, then that means z' = cos(100x^2 + z^3) > 0 (using the differential equation),
    # i.e. 100x^2 + z^3 is in an interval [-pi/2 + 2pik, pi/2 + 2pik] for some integer k.
    # If z^3 is decreasing, that means 100x^2 + z^3 is in an interval [pi/2 + 2pik, 3pi/2 + 2pik] for some integer k.
    # The lengths of these intervals are the same, but because 100x^2 >> z^3 in our range, the z^3 term can be thought of as "slightly speeding up or slowing down the 100x^2".
    # If the z^3 term is slowing things down, then 100x^2 + z^3 will linger in the interval above for a big longer,
    # whereas if the z^3 term is speeding things up, then it'll move past that interval more quickly and get to the next one over (of the other type) more quickly.
    # Thus, qualitatively we might expect that overall there will be this downward trend, since the function "spends more time in the decreasing zones than it does in the increasing zones".

    ### c)
    # The point of this question is to get you to push your computer as hard as you can, i.e. to get the most out of it.
    # My expectation is that this will manifest as things taking longer and longer, and at some point you just don't feel like waiting that long anymore.
    # In part a) I found that for deltax = 10^(-5), my code took about 1 second to run.
    # And when I assessed precision by taking deltax = 10^(-6), my code took about 10s to run.
    # I found that deltax = 10^(-5) was enough accuracy, and 10^(-3) wasn't.
    # Say I want to try to get xmax = 100.
    # Then 100x^2 = (100x)x will look like an angular frequency of order 10^4 around x = 100.
    # Taking deltax = 10^(-5) as above should be borderline for this.
    # Regarding runtime, because I'm taking ten times as many x values if I'm going up to xmax = 100 instead of xmax = 10,
    # I expect about 10s for deltax = 10^(-5) and 100s for deltax = 10^(-6).
    # Another problem is that without care, this kind of stuff can make your computer run out of ram:
    # If we store every value of z(x), that's 100/10^(-6) = 10^8 values, which has a memory requirement conveniently measured in gigabytes of RAM.
    # We can instead store every 100th value or similar, as well as the last one, and then we'll be fine.
    # I added the spacing optional variable to my Euler's method code get_soln to do this.
    # Let's try deltax = 10^(-5), and 10^(-6), because we expect these to run pretty quickly.
    xmax = 100
    print('For xmax = %i:' % xmax)
    start_time = time.time()
    z_vals_100_1e5 = get_soln(z0, x0, xmax, 10.0**(-5), F=F, spacing=1)
    runtime_1e5 = time.time() - start_time
    z100_val_1e5 = sorted(list(z_vals_100_1e5.items()))[-20:]
    print('last (x,z)s 1e5:')
    pprint.pprint(z100_val_1e5)
    print('runtime_1e5: ' + str(runtime_1e5) + ' seconds')
    # last (x,z)s 1e5:
    # [(99.99981002054179, 0.9293977351336113),
    #  (99.99982002054179, 0.9293879606189035),
    #  (99.9998300205418, 0.929377961436487),
    #  (99.9998400205418, 0.929368136119069),
    #  (99.9998500205418, 0.9293588762703148),
    #  (99.9998600205418, 0.9293505509645229),
    #  (99.9998700205418, 0.9293434920422701),
    #  (99.99988002054181, 0.9293379808867461),
    #  (99.99989002054181, 0.9293342372070649),
    #  (99.99990002054182, 0.9293324102755324),
    #  (99.99991002054182, 0.9293325729687908),
    #  (99.99992002054182, 0.9293347188516399),
    #  (99.99993002054183, 0.9293387624214821),
    #  (99.99994002054183, 0.9293445425055749),
    #  (99.99995002054183, 0.9293518286774434),
    #  (99.99996002054183, 0.9293603304380773),
    #  (99.99997002054184, 0.9293697087967041),
    #  (99.99998002054184, 0.9293795897895291),
    #  (99.99999002054184, 0.9293895793968073),
    #  (100.00000002054185, 0.9293992792622501)]
    # runtime_1e5: 10.337748527526855 seconds

    start_time = time.time()
    z_vals_100_1e6 = get_soln(z0, x0, xmax, 10.0**(-6), F=F, spacing=10)
    runtime_1e6 = time.time() - start_time
    z100_val_1e6 = sorted(list(z_vals_100_1e6.items()))[-20:]
    print('last (x,z)s 1e6:')
    pprint.pprint(z100_val_1e6)
    print('runtime_1e6: ' + str(runtime_1e6) + ' seconds')
    # last (x,z)s 1e6:
    # [(99.9998108395414, 0.9293801831870288),
    #  (99.99982083954137, 0.9293702561741799),
    #  (99.99983083954135, 0.9293603163016974),
    #  (99.99984083954132, 0.9293507597389917),
    #  (99.9998508395413, 0.9293419673854001),
    #  (99.99986083954127, 0.9293342896950135),
    #  (99.99987083954125, 0.9293280327119091),
    #  (99.99988083954122, 0.9293234458714574),
    #  (99.9998908395412, 0.9293207120536644),
    #  (99.99990083954117, 0.9293199402854776),
    #  (99.99991083954114, 0.9293211613840754),
    #  (99.99992083954112, 0.9293243267164181),
    #  (99.9999308395411, 0.929329310126396),
    #  (99.99994083954107, 0.9293359129545752),
    #  (99.99995083954104, 0.9293438719519562),
    #  (99.99996083954102, 0.929352869773217),
    #  (99.99997083954099, 0.9293625476313623),
    #  (99.99998083954097, 0.9293725196087965),
    #  (99.99999083954094, 0.9293823880531069),
    #  (100.00000083954092, 0.929391759442127)]
    # runtime_1e6: 71.42251706123352 seconds
    
    # Indeed we see that actually the values at 100 are pretty close.
    
    # Can we do xmax = 1000?
    # That probably requires deltax to be 10x smaller again (because the heuristic "pseudo angular frequency" thing will be 10x bigger at the end),
    # and it also requires 10x as many steps, because we're going 10 time farther.
    # That'll give a runtime 100x longer, i.e. 1000s if we take the deltax = 10^(-5) example above as a benchmark.
    # This is about 20 minutes, and I can easily run this in the background for that long while doing something else.
    # We shouldn't be too optimistic, because xmax = 100 with deltax = 10^(-5) was borderline.

    if not skip1k:
        xmax = 1000
        print('For xmax = %i:' % xmax)
        start_time = time.time()
        z_vals_1000_1e6 = get_soln(z0, x0, xmax, 10.0**(-6), F=F, spacing=100)
        runtime_1e6 = time.time() - start_time
        z1000_val_1e6 = sorted(list(z_vals_1000_1e6.items()))[-50:]
        print('last (x,z)s 1e6:')
        pprint.pprint(z1000_val_1e6)
        print('runtime_1e6: ' + str(runtime_1e6) + ' seconds')
        # For xmax = 1000:
        # last (x,z)s 1e6:
        # [(999.9951985672715, 0.9146974202038994),
        #  (999.9952985672712, 0.9146948739346374),
        #  (999.995398567271, 0.9146982323907404),
        #  (999.9954985672707, 0.9147035195776644),
        #  (999.9955985672705, 0.9147044764107528),
        #  (999.9956985672702, 0.9146999696094089),
        #  (999.99579856727, 0.9146953330022598),
        #  (999.9958985672697, 0.9146960540871532),
        #  (999.9959985672695, 0.914701278667042),
        #  (999.9960985672692, 0.9147048216868235),
        #  (999.996198567269, 0.9147024885531372),
        #  (999.9962985672687, 0.9146970402714554),
        #  (999.9963985672684, 0.9146949251643526),
        #  (999.9964985672682, 0.9146986462048913),
        #  (999.9965985672679, 0.9147037980585011),
        #  (999.9966985672677, 0.9147042816269907),
        #  (999.9967985672674, 0.9146995238046073),
        #  (999.9968985672672, 0.9146951557300084),
        #  (999.9969985672669, 0.9146963472308476),
        #  (999.9970985672667, 0.9147016872351951),
        #  (999.9971985672664, 0.9147048538393223),
        #  (999.9972985672662, 0.9147020978923034),
        #  (999.9973985672659, 0.9146966809849048),
        #  (999.9974985672657, 0.9146950145214813),
        #  (999.9975985672654, 0.914699070489892),
        #  (999.9976985672652, 0.9147040469245344),
        #  (999.9977985672649, 0.9147040521774444),
        #  (999.9978985672647, 0.9146990793332203),
        #  (999.9979985672644, 0.9146950142474685),
        #  (999.9980985672642, 0.9146966682686092),
        #  (999.9981985672639, 0.9147020827567669),
        #  (999.9982985672636, 0.9147048474161519),
        #  (999.9983985672634, 0.9147016887890381),
        #  (999.9984985672631, 0.9146963452510589),
        #  (999.9985985672629, 0.9146951415748867),
        #  (999.9986985672626, 0.91469950198272),
        #  (999.9987985672624, 0.9147042639425536),
        #  (999.9988985672621, 0.9147037895064082),
        #  (999.9989985672619, 0.9146986396116503),
        #  (999.9990985672616, 0.9146949098989849),
        #  (999.9991985672614, 0.9146970148772439),
        #  (999.9992985672611, 0.9147024619880665),
        #  (999.9993985672609, 0.9147048020950188),
        #  (999.9994985672606, 0.9147012642282631),
        #  (999.9995985672604, 0.9146960358310258),
        #  (999.9996985672601, 0.9146953055903698),
        #  (999.9997985672599, 0.9146999373181232),
        #  (999.9998985672596, 0.9147044470994342),
        #  (999.9999985672594, 0.9147034953392281),
        #  (1000.0000005672593, 0.9147018758647304)]
        # runtime_1e6: 679.1783785820007 seconds

        plot_1k = pts_plot_single(z_vals_1000_1e6)
        plot_1k_trunc = pts_plot_single(z_vals_1000_1e6, kmin=999)
        ga = graphics_array([plot_1k, plot_1k_trunc])
        ga.show(figsize=[16,10]) # a4f2.png
        # The plot looks funny here in part because because we're recording only every 100th value.
        # We'll see below that we get a similar-looking plot for deltax 10x smaller and recording every 1000th value

        # We should really do 10x, or at least 3x, the precision, i.e. deltax = 10^(-7).
        # We can still do this pretty easily by letting the code run in the background for 2-3 hours or so.
        # This serves as a check of the quality of the result above, but can also be used to give a sense of the accuracy of this answer.
        start_time = time.time()
        z_vals_1000_1e7 = get_soln(z0, x0, xmax, 10.0**(-7), F=F, spacing=1000)
        runtime_1e7 = time.time() - start_time
        z1000_val_1e7 = sorted(list(z_vals_1000_1e7.items()))[-1000:]
        print('last (x,z)s 1e7:')
        pprint.pprint(z1000_val_1e7)
        print('runtime_1e7: ' + str(runtime_1e7) + ' seconds')
        # For xmax = 1000:
        # last (x,z)s 1e7:
        # [(999.9001791478381, 0.9146946034841709),
        #  (999.9002791478038, 0.9146934857332912),
        #  (999.9003791477694, 0.9146881780361424),
        #  (999.900479147735, 0.9146849436581332),
        #  (999.9005791477007, 0.9146875990679251),
        #  (999.9006791476663, 0.9146930096791851),
        #  (999.900779147632, 0.914694789106447),
        #  ...
        #  (999.9984791140637, 0.9146898876047316),
        #  (999.9985791140293, 0.9146940857343889),
        #  (999.998679113995, 0.9146926406721593),
        #  (999.9987791139606, 0.9146872623364575),
        #  (999.9988791139263, 0.9146843166726419),
        #  (999.9989791138919, 0.9146872899982058),
        #  (999.9990791138575, 0.9146926618355368),
        #  (999.9991791138232, 0.914694072264266),
        #  (999.9992791137888, 0.9146898509229019),
        #  (999.9993791137545, 0.914684994248379),
        #  (999.9994791137201, 0.9146852507949221),
        #  (999.9995791136857, 0.9146903161994),
        #  (999.9996791136514, 0.9146941932241918),
        #  (999.999779113617, 0.9146922914518425),
        #  (999.9998791135827, 0.9146868615125846),
        #  (999.9999791135483, 0.9146843306562392),
        #  (1000.0000000135411, 0.9146914435352426)]
        # runtime_1e7: 7230.173061132431 seconds

        plot_1k = pts_plot_single(z_vals_1000_1e7)
        plot_1k_trunc = pts_plot_single(z_vals_1000_1e7, kmin=999)
        ga = graphics_array([plot_1k, plot_1k_trunc])
        ga.show(figsize=[16,10]) # a4f3.png

        # You can see that the z values above don't overlap with the ones for 10^(-6), 0.91469-0.91468 vs. 0.91470-0.91469, which is bad news.
        # Another thing you can notice in the above is that the errors from float addition are starting to add up.
        # We took deltax = 10^(-7), and then, over the course of the computation, added up deltax + deltax + ... + deltax 10^10 times.
        # You can see that all the x values near the end end in a bunch of random digits that are not 0.
        # This is another thing to worry about when taking this computation out this far.
        # My biggest concern here is that it is pretty hard to get a good sense of the accuracy of this result.
        # Getting a reasonable sense of the accuracy and precision of any experiment is very often more important than getting accuracy as high as possible.
        # For that reason, I wouldn't put too much stock into this xmax = 10^3 case.

        # Going to xmax = 10^4 probably requires another two orders of magnitude in the runtime.
        # Running it in the background for this long is realistic, but I would have had to start this assignment in advance.

    ### d)
    # Like in the previous problems, writing cos(100x^2 + z^3) = cos((100x)x + z^3) lets us think about y' as oscillating with angular frequency roughly 100x.
    # This can be made precise by comparing z'(x)*deltax = F(z,x) to 1/2*z''(x)*deltax^2 = d/dx F(z,x) = partial F/partial x + partial F/partial z * z'(x),
    # also as outlined in one of the previous solutions. This leads us to want to take deltax << 1/xmax.
    # 1. It's kind of tedious to have to pick deltax differently each time, test, etc.
    # Wouldn't it be nice if you could just have the computer pick the right size for you?
    # 2. Around, say, x = 100, you need 100x more precision, i.e. deltax 100x smaller, than around x = 1.
    # It's a bit of a waste to take the really small deltax for the whole thing.
    # So let's try scaling deltax like 1/x:
    def get_soln_dynamic(y0, t0, tmax, deltat_coeff, F=None, spacing=1, max_counter=10**10):
        y_vals = {t0:y0}
        y = float(y0)
        t = float(t0)
        counter = 0
        while (t < tmax) and (counter < max_counter):
            if F is None:
                dydt = get_dydt(y,t)
            else:
                dydt = F(y,t)
            deltat = deltat_coeff / max(1,t)
            new_t = t + deltat
            new_y = y + dydt * deltat
            t = new_t
            y = new_y
            if (counter%spacing) == 0:
                y_vals[t] = y
            counter += 1
        y_vals[t] = y
        return y_vals
    
    xmax = 10
    deltax_coeff = 10.0**(-4) # So that it's around 10^(-5) when x ~= 10, which we found in part a) was reasonable
    print('For xmax = %i:' % xmax + ' and dynamic step sizes with deltax_coeff = ' + str(deltax_coeff))
    start_time = time.time()
    z_vals_10_dynamic = get_soln_dynamic(z0, x0, xmax, deltax_coeff, F=F)
    runtime_10_dynamic = time.time() - start_time
    to_print = sorted(list(z_vals_10_dynamic.items()))[-4000:-20:100] + sorted(list(z_vals_10_dynamic.items()))[-20:]
    print('last (x,z)s xmax = %i dynamic:' % xmax)
    pprint.pprint(to_print)
    print('runtime_10_dynamic: ' + str(runtime_10_dynamic) + ' seconds')
    # For xmax = 10: and dynamic step sizes with deltax_coeff = 0.0001
    # last (x,z)s xmax = 10 dynamic:
    # [(9.959931217608144, 0.9444308541176435),
    #  (9.960935190511218, 0.9450730360299813),
    #  (9.961939062232938, 0.9443066456366978),
    #  (9.962942832803877, 0.9443003872882635),
    #  (9.963946502254597, 0.9450700677617414),
    #  (9.96495007061566, 0.9444337837307399),
    #  (9.965953537917612, 0.9441918528275709),
    #  (9.966956904190978, 0.9450275585216219),
    #  (9.967960169466245, 0.9445719114017853),
    #  (9.968963333773933, 0.9441137223640518),
    #  (9.969966397144498, 0.9449487382788648),
    #  (9.970969359608407, 0.944709882176684),
    #  (9.97197222119612, 0.9440720291815236),
    #  (9.97297498193807, 0.9448397362069472),
    #  (9.973977641864645, 0.9448365687063024),
    #  (9.974980201006266, 0.9440698922110671),
    #  (9.97598265939332, 0.944709093586473),
    #  (9.976985017056185, 0.9449417423806931),
    #  (9.977987274025198, 0.9441072752219212),
    #  (9.978989430330707, 0.9445670811075846),
    #  (9.979991486003035, 0.945016882930602),
    #  (9.980993441072497, 0.9441809984249527),
    #  (9.981995295569394, 0.9444248753878997),
    #  (9.982997049523982, 0.9450558549626791),
    #  (9.983998702966538, 0.944284999288179),
    #  (9.985000255927318, 0.9442936617923658),
    #  (9.986001708436548, 0.9450553979335141),
    #  (9.987003060524442, 0.9444108198767688),
    #  (9.988004312221214, 0.9441837371114503),
    #  (9.989005463557051, 0.9450153897823536),
    #  (9.990006514562118, 0.9445482809067421),
    #  (9.991007465266573, 0.9441036856112346),
    #  (9.99200831570058, 0.944938861583994),
    #  (9.993009065894244, 0.9446862890219762),
    #  (9.994009715877683, 0.9440596954253914),
    #  (9.995010265680994, 0.9448317600709263),
    #  (9.996010715334263, 0.9448137142958252),
    #  (9.997011064867554, 0.9440550699511255),
    #  (9.998011314310922, 0.9447024752694817),
    #  (9.999011463694414, 0.944920270048535),
    #  (9.999811511178256, 0.9442030650132973),
    #  (9.999821511366749, 0.9441959233621469),
    #  (9.99983151154524, 0.9441889230014537),
    #  (9.999841511713731, 0.9441820667283938),
    #  (9.999851511872222, 0.9441753572828243),
    #  (9.999861512020713, 0.9441687973461875),
    #  (9.999871512159203, 0.9441623895404374),
    #  (9.999881512287692, 0.944156136426989),
    #  (9.999891512406181, 0.9441500405056921),
    #  (9.99990151251467, 0.9441441042138285),
    #  (9.999911512613158, 0.9441383299251347),
    #  (9.999921512701647, 0.944132719948848),
    #  (9.999931512780135, 0.9441272765287794),
    #  (9.999941512848622, 0.9441220018424118),
    #  (9.99995151290711, 0.9441168980000237),
    #  (9.999961512955597, 0.9441119670438398),
    #  (9.999971512994085, 0.9441072109472082),
    #  (9.999981513022572, 0.9441026316138048),
    #  (9.999991513041058, 0.9440982308768653),
    #  (10.000001513049545, 0.9440940104984442)]
    # runtime_10_dynamic: 0.4980177879333496 seconds

    plot_dynamic = pts_plot_single(z_vals_10_dynamic, color='darkgoldenrod')
    plot_dynamic_zoomed = pts_plot_single({x:z for x,z in z_vals_10_dynamic.items() if x > 9.9}, color='darkgoldenrod')
    ga = graphics_array([plot_good, plot_bad, plot_dynamic, plot_good_zoomed, plot_bad_zoomed, plot_dynamic_zoomed], 2, 3)
    ga.show(figsize=[33,16]) # a4f11.png

    # Now we can change xmax and just run again, without having to think about deltax:
    xmax = 100
    print('For xmax = %i:' % xmax + ' and dynamic step sizes with deltax_coeff = ' + str(deltax_coeff))
    start_time = time.time()
    z_vals_100_dynamic = get_soln_dynamic(z0, x0, xmax, deltax_coeff, F=F)
    runtime_100_dynamic = time.time() - start_time
    to_print = sorted(list(z_vals_100_dynamic.items()))[-40000:-20:1000] + sorted(list(z_vals_100_dynamic.items()))[-20:]
    print('last (x,z)s xmax = %i dynamic:' % xmax)
    pprint.pprint(to_print)
    print('runtime_100_dynamic: ' + str(runtime_100_dynamic) + ' seconds')
    # For xmax = 100: and dynamic step sizes with deltax_coeff = 0.0001
    # last (x,z)s xmax = 100 dynamic:
    # [(99.95999329974262, 0.9293425668910429),
    #  (99.9609936949688, 0.9293719536301437),
    #  (99.96199408018326, 0.929425692439999),
    #  (99.96299445538631, 0.9294400827048727),
    #  (99.96399482057825, 0.9293980183590124),
    #  (99.96499517575937, 0.9293492194115683),
    #  (99.96599552093, 0.9293513747217148),
    #  (99.9669958560904, 0.9294018613666394),
    #  (99.9679961812409, 0.9294408297693262),
    #  (99.96899649638182, 0.9294220714476237),
    #  (99.96999680151342, 0.9293677222713188),
    #  (99.970997096636, 0.9293420397802742),
    #  (99.97199738174987, 0.9293753537030643),
    #  (99.97299765685534, 0.9294281499450756),
    #  (99.97399792195272, 0.9294378442594722),
    #  (99.97499817704224, 0.9293928903481466),
    #  (99.97599842212429, 0.9293464278542662),
    #  (99.97699865719913, 0.9293533805476629),
    #  (99.97799888226709, 0.9294054460358135),
    #  (99.97899909732843, 0.9294409050441155),
    #  (99.97999930238348, 0.9294177048374805),
    #  (99.98099949743252, 0.9293632389388481),
    #  (99.98199968247587, 0.9293419029301018),
    #  (99.98299985751379, 0.9293788818793278),
    #  (99.98400002254662, 0.9294303218684119),
    #  (99.98500017757468, 0.9294352449296801),
    #  (99.98600032259819, 0.9293877534854554),
    #  (99.98700045761751, 0.9293439901302302),
    #  (99.98800058263292, 0.9293556842806723),
    #  (99.98900069764476, 0.9294089200970224),
    #  (99.9900008026533, 0.929440592389311),
    #  (99.99100089765882, 0.9294131324388919),
    #  (99.99200098266164, 0.929358975661636),
    #  (99.99300105766206, 0.9293421517513278),
    #  (99.99400112266038, 0.9293825049234609),
    #  (99.99500117765689, 0.9294321856896394),
    #  (99.99600122265193, 0.9294322995674894),
    #  (99.99700125764575, 0.9293826423999698),
    #  (99.99800128263865, 0.9293419196696762),
    #  (99.99900129763095, 0.9293582622402186),
    #  (99.99998130262118, 0.9293930579302179),
    #  (99.99998230262136, 0.9293940553916669),
    #  (99.99998330262154, 0.9293950512292561),
    #  (99.99998430262171, 0.9293960450445614),
    #  (99.99998530262187, 0.929397036439968),
    #  (99.99998630262202, 0.9293980250188297),
    #  (99.99998730262216, 0.929399010385628),
    #  (99.99998830262228, 0.9293999921461299),
    #  (99.9999893026224, 0.9294009699075461),
    #  (99.9999903026225, 0.929401943278688),
    #  (99.9999913026226, 0.9294029118701245),
    #  (99.99999230262269, 0.9294038752943374),
    #  (99.99999330262277, 0.9294048331658773),
    #  (99.99999430262284, 0.9294057851015171),
    #  (99.99999530262289, 0.9294067307204058),
    #  (99.99999630262293, 0.9294076696442211),
    #  (99.99999730262297, 0.9294086014973201),
    #  (99.999998302623, 0.9294095259068906),
    #  (99.99999930262301, 0.9294104425030997),
    #  (100.00000030262302, 0.9294113509192419)]
    # runtime_100_dynamic: 66.9851906299591 seconds

    plot_100_dynamic = pts_plot_single(z_vals_100_dynamic, plot_spacing=10, color='darkgoldenrod')
    plot_100_dynamic_zoomed = pts_plot_single(z_vals_100_dynamic, kmin=99.99, color='darkgoldenrod')
    ga = graphics_array([plot_100_dynamic, plot_100_dynamic_zoomed])
    ga.show(figsize=[16,10]) # a4f12.png # This takes a while

    # A dynamic step size emphasizes how overall these numerical calculations take longer and longer as xmax increases:
    # The number of steps to get from x = n to x = n+1 is approximately n / deltax_coeff.
    # That means the number of steps to get all the way to xmax is roughly xmax^2 / 2deltax_coeff.
    # This factor of 2 in the denominator means that we should expect that this will be about twice as fast as constant step size.
    # We saw that we did save a factor of 2 for xmax=10: 1s vs. 0.5s.
    # For xmax = 100 the times were similar though. Probably because of the additional overhead of maintaining the dynamic stepsizes.
    # Moreover, we got slightly different answers in each case.
    # My takeaway is that overall dynamic stepsize is a bit more convenient here, but doesn't make much of a concrete difference.

    # The point of asking about cos(100x^2.1 + z^3) is the following.
    # One thing someone might try for dynamic stepsizes is to pick a step size as a function of how many steps have been taken so far.
    # This is a bad idea if done poorly.
    # If one takes deltax_n = constant/n on the nth step, then, while 1/1 + 1/2 + 1/3 + ... + 1/n ~= log(n) -> infinity,
    # it goes to infinity very slowly: it'll take exponentially many steps to reach xmax.
    # For cos(100x^2.1 + z^3), if one tries deltax_n = constant/n^1.1, then the sum of step sizes converges, i.e. there are certain xmax's you'd never reach.
    # You can mirror the stepsize constant/x for the x^2 case that we did above by picking deltax_n = constant/sqrt(n), this is a little sneaky and overcomplicated.

    return



def problem_2():
    ### a)
    # I'll do this using complex exponentials. I explain how to switch to cos at the end.
    # Given the equation y'' + 2cy' + d^2y = f*exp(ikt), we can guess y = Aexp(ikt) as a particular solution.
    # Substituting gives (-k^2 + i2ck + d^2)Aexp(ikt) = f*exp(ikt)
    # and solving for A gives A = A(k) = f/(d^2 - k^2 + i2ck)
    # To find the resonance frequency, we want to find k such that |A(k)| is maximized.
    # It's equivalent, and less work, to minimize |A(k)|^(-2) = (d^2 - k^2)^2 + 4c^2k^2
    # Differentiating with respect to k gives d/dk |A(k)|^(-2) = -4k(d^2 - k^2) + 8kc^2
    # To find the minimum, we set this derivative to 0 and solve for k.
    # The value k = 0 is always a zero of the derivative, and not the one we're looking for. Let's assume k > 0 so we can divide by 4k.
    # We get 2c^2 + k^2 - d^2 = 0, i.e. k = sqrt(d^2 - 2c^2).
    # This is what we got in class. We had different variable names: k = omega, d = omega_0, c = lambda.
    # For 2c = 0.1 and d^2 = 2, we find that the resonance frequency is sqrt(2 - 0.1^2/2) = sqrt(1.995) ~= 1.412444689...
    # The corresponding amplitude is A(sqrt(d^2 - 2c^2)) = f/(d^2 - sqrt(d^2 - 2c^2)^2 + i2c*sqrt(d^2 - 2c^2)) = f/(2c^2 + i2c*sqrt(d^2 - 2c^2)).
    # In this case, f = 1, and A(sqrt(1.995)) = 1/(0.005 + i*0.1*sqrt(1.995)).
    # The amplitude is the absolute value of A, which here is 1/|0.005 + i*0.1*sqrt(1.995)|
    # = 1/sqrt(0.005^2 + 0.01*1.995) = 1/sqrt(0.019975) ~= 7.07549...
    #
    # To turn the right hand side into cos, let's take real parts: Re(y'' + 2cy' + d^2y) = Re(f*exp(ikt))
    # Because the coefficients in the differential equation are real, we have Re(y'') + 2c*Re(y') + d^2*Re(y) = f*cos(kt)
    # Writing y(t) = u(t) + iv(t) for u,v real functions, we see Re(y') = Re(u' + iv') = u' = Re(y)'
    # So the real part of y solves the differential equation with Re(f*exp(ikt)) = f*cos(kt) on the right.
    # This is y(t) = Re(A(k) * exp(ikt)) = Re(f(cos(kt) + isin(kt))/(d^2 - k^2 + i2ck))
    # = Re(f/((d^2 - k^2)^2 + 4c^2k^2) * (cos(kt) + isin(kt)) * (d^2 - k^2 - i2ck))
    # = f*((d^2 - k^2)*cos(kt) + 2ck*sin(kt))/((d^2 - k^2)^2 + 4c^2k^2).
    # We can then plug in the values d^2 = 2, k = sqrt(1.995), c = 0.05, and f = 1 to get the particular solution.
    # As explained in Wainright-West around equation (2.85), the amplitude with just cos on the right hand side is the same.

    ### b)
    # This question is intended to highlight that "initial conditions such that the homogeneous part is 0" is not the same as "y(o) = y'(o) = 0".
    # In the former case, from the discussion above, we have
    # y(t) = f*((d^2 - k^2)*cos(kt) + 2ck*sin(kt))/((d^2 - k^2)^2 + 4c^2k^2)
    # = (0.005*cos(sqrt(1.995)t) + 0.1sqrt(1.995)*sin(sqrt(1.995)) / (0.005^2 + 0.01*1.995)
    # This, like all linear combinations of sin and cos with the same frequencies, is a pure sine wave, and reaches its maximum amplitude immediately.
    # I suspect many students might doubt themselves answering this, and it's very important to be confident.
    
    # The initial conditions that the given ones contrast with are y(0) = y'(0) = 0. Let's see what happens in that case.
    def d2ydt2(y, yp, t, omega):
        to_ret = -0.1*yp - 2*y + cos(omega*t)
        return to_ret
    
    def get_soln_omega(y0, yp0, t0, tmax, deltat, omega):
        y_vals = {t0:y0}
        y = y0
        yp = yp0
        t = t0
        while t < tmax:
            new_yp = yp + d2ydt2(y, yp, t, omega) * deltat
            new_y = y + yp * deltat
            new_t = t + deltat
            yp = new_yp
            y = new_y
            t = new_t
            y_vals[t] = y
        return y_vals

    def pts_plot_omega_range(omegamin, omegamax, num_omegas, y0, yp0, t0, tmax, deltat, spacing=1, kmin=None, kmax=None, vmin=None, vmax=None, absvmax=None, **kwargs):
        deltaomega = float(omegamax - omegamin)/num_omegas
        omegalist = [omegamin + deltaomega*i for i in range(num_omegas+1)]
        soln_dict = {omega:get_soln_omega(y0, yp0, t0, tmax, deltat, omega) for omega in omegalist}
        pts_list = [pts_plot_single(soln_dict[omega], spacing=spacing, kmin=kmin, kmax=kmax, vmin=vmin, vmax=vmax, absvmax=absvmax, color=hue(0.32+0.9*i/len(soln_dict))) for i,omega in enumerate(sorted(list(soln_dict.keys())))]
        return pts_list

    tmax = 100
    pts_list = pts_plot_omega_range(1, 3, 100, 0, 0, 0, tmax, 0.001)
    show(sum(pts_list)) # a4f4.png

    cutoff = 0.99*1/sqrt(0.005**2 + 0.01*1.995)
    resonance_soln = get_soln_omega(0, 0, 0, tmax, 0.001, sqrt(1.995))
    pts_resonance = pts_plot_single(resonance_soln)
    t99 = sorted(list([k for k,v in resonance_soln.items() if abs(v) > cutoff]))[0]
    cutoff_line_plus = line([(0, cutoff), (tmax, cutoff)], color='black', dpi=400)
    cutoff_line_minus = line([(0, -cutoff), (tmax, -cutoff)], color='black', dpi=400)
    pt_99 = points([(t99, cutoff)], size=10, color='red', faceted=True, markeredgecolor='black', zorder=99)
    show(pts_resonance + cutoff_line_plus + cutoff_line_minus + pt_99) # a4f5.png
    print(t99) # 72.23599999998086

    # If you imagine yourself on a playground swing set, starting at rest and swinging under your own power,
    # then that's reasonably modelled by something like the current ODE with y(0) = y'(0) = 0.
    # You know intuitively that you ramp up to max amplitude over many swings.
    # Because we said above that the particular solution with homogeneous part 0 is a pure sine wave,
    # that means that the ramping up on a swing set is in fact the homogeneous solution decaying!
    # Personally I find this quite weird and cool.
    
    return



def problem_3():
    def F(y,t):
        try:
            to_ret = -y*t + 0.1*y**3
        except OverflowError: # y' = y grows exponentially, so y' = y^3 grows really really fast. We don't care exactly how big it gets.
            to_ret = 10**6
        return to_ret

    # Let's start by just looking at the solutions, and assess from there
    plt_list = pts_plot_range(0, 4, 1000, 0, 6, 0.001, F=F, vmax=10)
    show(sum(plt_list)) # a4f6.png
    # We see that the critical alpha will be between 2 and 3, and the converging solutions don't go above 5.
    
    def get_critical_alpha(alpha_min, alpha_max, eps, deltat):
        # binary search
        counter = 0
        max_counter = 100
        while ((alpha_max - alpha_min) > eps) and (counter < max_counter):
            counter += 1
            alpha = (alpha_min + alpha_max)/2.0
            soln = get_soln(alpha, 0, 4, deltat, F=F)
            if any([v > 10 for v in soln.values()]):
                alpha_max = alpha
            else:
                alpha_min = alpha
        if counter >= max_counter:
            raise ValueError
        return (alpha_min + alpha_max)/2.0
    
    alpha0 = get_critical_alpha(2, 3, 10.0**(-6), 10.0**(-6)) # ~1 minute
    print(alpha0) # 2.375267505645752
    print(sqrt(10.0)/float(pi)**0.25) # 2.37526752924330
    
    return



def problem_4():
    def F(y,t):
        return 3*t**2 / (3*y**2 - 4.0)
    # Let's start by just looking at the solutions
    t0 = 0
    tmax = 6
    deltat = 0.01
    plt_list = pts_plot_range(-5, 5, 1000, t0, tmax, deltat, F=F) # this is 10^3 y0 values and about 10^3 steps for each, so this should take a couple seconds to run
    show(sum(plt_list)) # a4f7.png
    # We see that we get something totally bizarre.
    # Let's cap the range of y values being plotted and try again:
    plt_list = pts_plot_range(-5, 5, 1000, t0, tmax, deltat, F=F, absvmax=10)
    show(sum(plt_list)) # a4f8.png
    # By the mixing of the colours, we can see that a lot of the solutions are jumping abruptly.
    # With that in mind, the previous image makes more sense: those random horizontal lines were solutions that jumped really far.
    # We see some obvious divergence for y a bit larger than 1,
    # and some weird noise for y a bit smaller than -1.

    # What's causing the solutions to jump?
    # Looking at the differential equation, we see that when y = +-sqrt(4/3), there will be division by 0.
    line_plus = line([(0,sqrt(4.0/3)), (tmax,sqrt(4.0/3))], dpi=400, color='black')
    line_minus = line([(0,-sqrt(4.0/3)), (tmax,-sqrt(4.0/3))], dpi=400, color='black')
    show(sum(plt_list) + line_plus + line_minus) # a4f9.png
    # For example, if you follow the dark blue solutions starting with y0 just above 0, they head downwards to the line -sqrt(4/3),
    # and when they reach that line, that's when you start seeing dark blue lines towards the top and bottom of the plot.
    # These are solutions that jumped because y' = F(y,t) was enormous at one of the steps, because of near division by 0.
    
    # At the upper line y = sqrt(4/3), things aren't so bad, because solutions diverge away from that line.
    # You can think of a direction field: above the line, F(y,t) > 0, so solutions keep increasing,
    # and just below y = sqrt(4/3) we have F(y,t) < 0, so the solutions decrease further.
    # The reason you see some solutions seemingly popping up from nowhere around this line is because solutions from near the bottom line might randomly jump there.
    # This is of course one of the many numerical artifacts, and doesn't actually happen.
    
    # The behaviour around the bottom line is much more insidious.
    # The line y = -sqrt(4/3) attracts nearby solutions,
    # because solutions just above the line will have negative derivative, and solutions just below the line will have positive derivative.
    # This is why, e.g. you see a blue solution at the very top of the plot. It started above the line, jumped just below it, and then got really close and jumped way up.
    # This is also why you see noise that looks like turbulance along the length of the line.

    # Changing the step size, we see that the details of the plot change completely:
    deltat = 0.001
    plt_list = pts_plot_range(-5, 5, 1000, t0, tmax, deltat, F=F, absvmax=10, spacing=10) # I'm doing spacing=10 here so that the plots are more easily comparable
    show(sum(plt_list) + line_plus + line_minus) # a4f10.png
    # You can see that all the spurrious curves are in different places and are different colours

    # This differential equation is separable:
    # int 3y^2 - 4 dy = int 3t^2 dt <=> y^3 - 4y = t^3 + C
    # It's possible to solve for y here, using the cubic formula, but very gross;
    # see https://www.wolframalpha.com/input?i=y%5E3+-+4y+%3D+x or https://www.wolframalpha.com/input?i=y%27+%3D+3t%5E2%2F%283y%5E2+-+4%29
    # You could solve for t, and get t = y(1 - 4/y^2 + C/y^3)^(1/3)
    # This shows that for y and t positive and large, the function is something like a line. This can be seen on the plot.
    # You can also go plot this function: https://www.wolframalpha.com/input?i=x%281+-+4%2Fx%5E2%29%5E%281%2F3%29&assumption=%22%5E%22+-%3E+%22Real%22
    # You see that for y around -sqrt(4/3) there is a bump, and t does not extend above that.
    
    return




#This is the standard boilerplate that calls the main() function.
if __name__ == '__main__':
    if '-profile' in sys.argv:
        cProfile.run('main()')
    else:
        main()

