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]: