diff --git a/src/core/time_arrays.jl b/src/core/time_arrays.jl
index da6e1890d..e25efc5f6 100644
--- a/src/core/time_arrays.jl
+++ b/src/core/time_arrays.jl
@@ -127,6 +127,59 @@ function _get_ts_indices(ts_array::Array{TimestepValue{T}, 1}, times::Union{Tupl
 	return [_get_time_value_position(times, ts) for ts in ts_array]
 end
 
+# Base.firstindex and Base.lastindex
+function Base.firstindex(arr::TimestepArray{T_TS, T, N, ti}) where {T_TS, T, N, ti}
+	if ti == 1
+		return Mimi.TimestepIndex(1)
+	else
+		return 1
+	end
+end
+
+function Base.lastindex(arr::TimestepArray{T_TS, T, N, ti}) where {T_TS, T, N, ti}
+	if ti == length(size(arr.data))
+		return Mimi.TimestepIndex(length(arr.data))
+	else
+		return length(arr.data)
+	end
+end
+
+function Base.lastindex(arr::TimestepArray{T_TS, T, N, ti}, dim::Int) where {T_TS, T, N, ti}
+	if ti == dim
+		return Mimi.TimestepIndex(size(arr.data, dim))
+	else
+		return size(arr.data, dim)
+	end
+end
+
+function Base.firstindex(arr::TimestepArray{T_TS, T, N, ti}, dim::Int) where {T_TS, T, N, ti}
+	if ti == dim
+		return Mimi.TimestepIndex(1)
+	else
+		return 1
+	end
+end
+
+# add axes methos copied from abstarctarray.jl:56
+function Base.axes(A::TimestepArray{T_TS, T, N, ti}, d::Int) where {T_TS, T, N, ti}
+	_d_lessthan_N = d <= N;
+	if d == ti
+		if _d_lessthan_N
+			return Tuple(TimestepIndex.(1:size(A,d)))
+		else
+			return TimestepIndex(1)
+		end
+	else
+		if _d_lessthan_N
+			if _d_lessthan_N
+				return 1:size(A,d)
+			else
+				return 1
+			end
+		end
+	end
+end
+
 #
 # b. TimestepVector
 #
@@ -233,8 +286,6 @@ function Base.length(v::TimestepVector)
 	return length(v.data)
 end
 
-Base.lastindex(v::TimestepVector) = length(v)
-
 #
 # c. TimestepMatrix
 #
@@ -267,7 +318,6 @@ function Base.getindex(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, 2}, id
 end
 
 function Base.getindex(mat::TimestepMatrix{VariableTimestep{TIMES}, T, 2}, idx::AnyIndex, ts::VariableTimestep{TIMES}) where {T, TIMES}
-	# WAS THIS: data = mat.data[ts.t, idx, ts.t]
 	data = mat.data[idx, ts.t]
 	_missing_data_check(data, ts.t)
 end
diff --git a/test/test_timesteparrays.jl b/test/test_timesteparrays.jl
index 8160d89c6..5a93a1aac 100644
--- a/test/test_timesteparrays.jl
+++ b/test/test_timesteparrays.jl
@@ -44,7 +44,7 @@ idx4 = TimestepIndex(4)
 
 x = TimestepVector{FixedTimestep{2000, 1}, Int}([9, 10, 11, 12])
 @test length(x) == 4
-@test lastindex(x) == 4
+@test lastindex(x) == TimestepIndex(4)
 
 time_dim_val = [9, 10, 11, 12]
 temp_dim_val = [100, 101, 102, 103]
@@ -611,6 +611,23 @@ y_mat = TimestepMatrix{VariableTimestep{y_years}, Int, 1}(time_dim_val[:,:,2])
 
 @test eltype(x_vec) == eltype(y_vec) == eltype(y_vec) == eltype(y_mat) == eltype(time_dim_val)
 
+# TODO begin syntax is depreacated v1.0.0 - v1.3.0, so enable the tests after 
+# after we disable Julia versions below v1.4.0
+
+# @test x_vec[begin] == time_dim_val[:,1,1][begin]
+# @test x_mat[begin,1] == time_dim_val[:,:,1][begin,1]
+# @test x_mat[begin,2] == time_dim_val[:,:,1][begin,2]
+# @test y_vec[begin] == time_dim_val[:,2,2][begin]
+# @test y_mat[begin,1] == time_dim_val[:,:,2][begin,1]
+# @test y_mat[begin,2] == time_dim_val[:,:,2][begin,2]
+
+@test x_vec[end] == time_dim_val[:,1,1][end]
+@test x_mat[end,1] == time_dim_val[:,:,1][end,1]
+@test x_mat[end,2] == time_dim_val[:,:,1][end,2]
+@test y_vec[end] == time_dim_val[:,2,2][end]
+@test y_mat[end,1] == time_dim_val[:,:,2][end,1]
+@test y_mat[end,2] == time_dim_val[:,:,2][end,2]
+
 #------------------------------------------------------------------------------
 # 6. Test that getindex for TimestepArrays doesn't allow access to `missing`
 #       values during `run` that haven't been computed yet.