In [51]:
using Plots; plotly(size=(200,200)); plot(Plots.fakedata(5,10))
;
In [52]:
# Define the ODE for Hooke's law
k = 0.5
s = [0; 1.0]
f = (t,y) -> [0 1.0; -k 0.0]*y

# Do simulation
h = 0.1
Tmax = 100
tsteps = collect(0.:h:Tmax)
nsteps = length(tsteps)
y = s # this is our current time
yhist = zeros(nsteps,2)
yhist[1,:] = y
for i=2:nsteps
    y = y + h*f(tsteps[i], y)
    yhist[i,:] = y
end
plot(tsteps, yhist[:,1])
Out[52]:
In [55]:
# Define the ODE for Malthusian Pop. Growth
# And use Heun's method
lam = -10.0
s = 10000.0
f = (t,y) -> lam*y

# Do simulation
#h = abs((2+0.01)/lam) # won't work
h = abs((2-0.01)/lam) # will work
# h = 0.1
Tmax = 100.
tsteps = collect(0.:h:Tmax)
nsteps = length(tsteps)
y = s # this is our current time
yhist = zeros(nsteps)
yhist[1] = y
for i=2:nsteps
    q1 = f(tsteps[i], y)
    q2 = f(tsteps[i]+h, y+h*q1) 
    y = y + h/2*(q1+q2)
    yhist[i] = y
end
plot(tsteps, yhist)
Out[55]:

Show the region of absolute stability for Heun's method

In [50]:
#using Images
nwidth = 300
nheight = 200
xs = linspace(-4,2,nwidth)
ys = linspace(-2,2,nheight)

X = zeros(nwidth,nheight)

for i=1:nwidth
    for j=1:nheight
        z = xs[i] + ys[j]*1im
        X[i,j] = abs(1+z+z*z/2.) < 1
    end
end

heatmap(xs,ys,X')
Out[50]: